有限元分析与数值计算：FEM不是计算器，从物理建模到数值方案有三次选择权

把有限元分析当计算器用的一个标志性习惯是：画完几何体，选默认单元，点”全部生成”，跑。结果出来了一组彩色云图，觉得”有结果了就对了”。

有限元分析本质上是在解一个近似问题。物理上的偏微分方程在连续域上有无穷多个解，有限元把无穷维的函数空间投影到一个有限维的子空间——用分片多项式去逼近真实解。这个”逼近”的质量取决于三次关键的选择。

abstract globe with sunrise in space in vector eps 10

第一次选择：单元类型——不止是六面体和四面体的区别

一个做压力容器应力评定的项目用的是Abaqus的C3D8R（8节点六面体缩减积分单元）。应力评定中有一条关键路径——接管与筒体连接处的应力分类线（SCL，Stress Classification Line）。沿SCL做应力线性化，分解出薄膜应力和弯曲应力。

C3D8R是线性单元，一个单元内的位移场是一次多项式。SCL穿过的单元壁厚方向只有一层，意味着弯曲应力在单元内的分布被假设为线性——但接管根部壁厚方向的实际弯曲应力分布有至少二阶的分量（壳体的边缘弯曲效应）。

换成C3D20R（20节点二次六面体），弯曲应力沿壁厚的分布从线性变成二次——SCL分解出的弯曲应力从187MPa变到212MPa，差了13%。给定Q345R在300°C的许用弯曲应力是约185MPa（按ASME Section VIII Div 2），187MPa刚好在许用范围内，212MPa直接超了。一个单元类型的差异可以颠倒合格与否的判定。

判断单元够不够级的不二标准是做了p型收敛——在同一个网格上把单元从线性升到二次看结果变不变。壁厚方向需要至少两个二次单元才能抓准边缘弯曲效应的三次应力分量。

第二次选择：积分方案——减缩积分的沙漏模式

C3D8R中的”R”是减缩积分（Reduced Integration）——线性六面体在完全积分下需要2×2×2=8个高斯积分点，减缩积分只用1个。省了计算量，但引入了一个隐患：沙漏模式（Hourglassing）。

沙漏是零能变形模式——单元可以发生一种变形（节点位移非零）但单一积分点处的应变为零。这种变形在物理上不应该存在，是减缩积分的数值假象。

这个项目在法兰螺栓孔附近用了C3D8R——孔周围的应力梯度和弯曲变形正好会激发沙漏模式。Abaqus默认提供了一个小的沙漏刚度（默认约0.05倍的弹性刚度），但这个刚度如果不调到足够大，沙漏能会增长。后处理里检查沙漏能与内能的比值——法兰孔周围几个单元的沙漏能/内能比到了8%。超过5%就不能接受。

切换成C3D8I（非协调模式单元，对弯曲有更好的响应，且不加沙漏控制）或直接用C3D20R，沙漏能从总内能的8%降到了1%以下。

第三次选择：收敛准则——力残差和位移残差不是同一张试卷

Abaqus默认的收敛准则是力残差的相对范数<0.5%且位移校正的相对范数<1%。这两个准则在大多数问题里等价，但在接触占主导的问题里常常打架。

这个项目的法兰预紧力步里，力残差在第三个牛顿迭代就收敛了（0.3%），但位移校正还差2.1%——螺栓杆在法兰孔里还在微滑，位移没稳定下来。如果只用力残差作为收敛判据（有些工程师会调高位移容差），得到的接触状态是一个还没稳定的中间状态。

多加了两步牛顿迭代，位移校正收敛到0.7%——螺栓杆的接触压强分布从不稳定的局部高压变成了一个更平滑的分布，峰值接触压力从420MPa降到了370MPa。直接用前三个迭代的结果，螺栓的接触应力被高估了13%。

有限元分析的”黑箱”其实透明。你做每一次默认选择的时候，程序的自动决策可能已经悄悄偏离了物理现实。三次选择——单元类型、积分方案、收敛准则——决定了你离工程实际有多近。