在自动剖分网格中,三角形单元是最常使用的单元。四面体单元是三维默认。几何形状越复杂,三角形单元就越多。
三角形单元与四边形单元不是相互等价的选项。它们在弯曲问题中的计算存在显著的数量级差异。
本文的对比原则是公平。两组单元使用完全相同的节点集。差异只来自单元本身。
三节点三角形单元的应变矩阵是常数,见式 (1)。
四节点四边形单元的应变随坐标线性变化。它可以用一个单元表示弯曲变形的线性应力梯度。三角形的梯度只能以马赛克形式拼接,见图 1。
错误选择的后果有两个:
(1) 弯曲刚度被高估。线性三角形比精确解偏硬。
(2) 应力收敛极慢。即使加密网格,单元中心的常应变值与真实应力之间的差缩小得极慢。这是式 (2) 中收敛阶的差异。
四边形单元在能量范数下
悬臂梁 L = 10,H = 1,端部加载抛物线剪力。同一节点集分别剖分为四边形单元与三角形单元,见图 2。
参考解为 Timoshenko梁的解析解:
误差定义为:
表 1 与图 3(a) 给出计算结果。
表 1 悬臂梁弯曲:同一节点集下的单元对比
从计算结果可以得到以下结论:
(1) 粗网格下三角形挠度偏低 46%,应力误差超过 100%。这不是网格质量问题。是常应变的先天限制。
(2) 三角形的应力误差几乎不降。节点数扩大 43 倍,应力误差仅从 108% 降到 22%。图 4 的跨中应力剖面给出原因。三角形单元只能给出阶梯状的常应变块,与真实线性分布始终存在结构性偏差。
(3) 达到 quad4 两个单元厚度的精度(约 −11%),三角形需要约 6 倍节点数。计算成本优势荡然无存。
在同一节点集下,四边形剖分为三角形的对角线方向会影响结果。nely = 4 网格下,统一对角线方向挠度误差为 −17.96%,交替翻转后 为−13.36%。对角线方向本身就是误差源的一部分。在使用自动剖分时应注意此现象。
考虑四分之一孔板。极坐标 O 网格 48×24,外边界均布拉应力,孔半径 0.5。四边形 1152 单元,三角形 2304 单元。结果见表 2。
表 2 孔板应力集中对比
有限宽修正后的理论峰值约 2.4 至 2.7。三角形的峰值落入理论区间,四边形因单元中心采样偏低。但三角形的两个对角线方向给出不同结果。剖分方向引入了额外的不确定性,见图 3(b)。
三角形在应力集中问题上并非不可用。但是需要考虑两个问题:
(1) 三角形单元的计算结果依赖剖分方向。
(2) 三角形单元需要用两倍的单元数才达到与四边形单元同等的网格密度。
在关心效率与可重复性的场合,四边形仍然占优。
综合两组实验的结果,单元选择准则见表 3。
表 3 三角形与四边形的选择准则
三维的对应关系明确。线性四面体单元与线性三角形单元共享同样的常应变缺陷。二阶四面体单元才有与线性六面体单元竞争的资格。商用软件的默认的二阶四面体选项正是为此准备。
线性三角形的缺陷源于常应变。其无法在一个单元内表示弯曲的应力梯度。使用后果是,弯曲刚度被高估与应力收敛停滞。在本文同一节点集的对比实验中,三角形单元应力误差是四边形单元的 10 到 100 倍。追求弯曲精度时应使用四边形单元或二阶单元。自动划分产生的三角形单元应被视作过渡手段,而非主力单元。