NAFEMS LE1 椭圆膜的内压问题 — V&V结果总结报告
V&V结果总结
老师，NAFEMS LE1的分析结果出来了，但我不太明白结果怎么看。我的参考解是92.7 MPa，但我的结果是92.5 MPa，这算通过吗？
好问题。LE1是椭圆膜受内压时的应力问题，评估点D处的 $\sigma_{yy}$ 参考解是 92.7 MPa。你的结果误差是 $(92.7 - 92.5)/92.7 \approx 0.22\%$，精度很高。但重点不只是单一结果，还要检查当网格细化时的收敛行为。
收敛行为是指改变网格的几个级别，看结果怎么变化吗？
完全正确。当网格逐级变细时，如果解单调收敛到参考解，说明结果的信任度高。反过来，如果数值在变化过程中波动，说明模型可能有问题。这次我们用7个求解器进行了验证，一起看看结果吧。
验证的背景和问题设置
LE1到底是什么问题？就是给椭圆膜加压力吗？
LE1是一个有椭圆形孔洞的板（平面应力状态）承受均匀内压的问题。椭圆的长轴和短轴比例会影响应力集中，评估点的 $\sigma_{yy}$ 有解析解，所以很适合用来检验FEA求解器的精度。比如说，汽车车身有椭圆形的开口（窗户或车门孔洞）时，能否正确计算周边的应力集中，这个问题就是测试这个能力。
全求解器结果的综合比较
用7个求解器计算了，结果怎么样？
用足够细的网格（48×32分割，二阶单元）计算的结果在这里。参考解是 $\sigma_{yy} = 92.70$ MPa。
| 求解器 | 单元类型 | 网格 | $\sigma_{yy}$ (MPa) | 与参考解的差 | 误差(%) | 判定 |
|---|---|---|---|---|---|---|
| MSC Nastran | CQUAD8 | 48×32 | 92.68 | -0.02 | 0.02 | 通过 |
| NX Nastran | CQUAD8 | 48×32 | 92.69 | -0.01 | 0.01 | 通过 |
| Abaqus | CPS8R | 48×32 | 92.71 | +0.01 | 0.01 | 通过 |
| Ansys | PLANE183 | 48×32 | 92.68 | -0.02 | 0.02 | 通过 |
| COMSOL | 二阶三角形 | 自动 | 92.72 | +0.02 | 0.02 | 通过 |
| Code_Aster | QUAD8 | 48×32 | 92.66 | -0.04 | 0.04 | 通过 |
| CalculiX | C2D8 | 48×32 | 92.63 | -0.07 | 0.08 | 通过 |
所有求解器的误差都不到0.1%？商用软件和开源软件几乎没有差别…
对，LE1是线性弹性的平面应力问题，求解器实现正确的话在理论上会收敛到相同的解。求解器间的偏差最多0.09 MPa（误差不足0.1%），这源于数值积分精度或单元定义的细微差异。实践中，「用哪个求解器」远不如「怎么划分网格」重要。从这些数据就能清楚看出来。
网格收敛行为的详细分析
网格粗细对结果的影响有多大？实践中往往没有时间把网格细化太多。
好问题。按单元类型用3个网格密度等级做对比，结果如下。
| 单元类型 | 粗网格误差(%) | 中网格误差(%) | 细网格误差(%) | 推荐 |
|---|---|---|---|---|
| QUAD4（一阶四边形） | 4.85 | 1.83 | 0.65 | 细网格及以上 |
| QUAD8（二阶四边形） | 0.43 | 0.11 | < 0.05 | 中网格及以上 |
| TRIA3（一阶三角形） | 15.2 | 6.69 | 2.10 | 不推荐 |
| TRIA6（二阶三角形） | 0.97 | 0.22 | 0.06 | 中网格及以上 |
TRIA3的粗网格误差15%，而QUAD8的粗网格只有0.43%！？差别这么大…
是的，这很令人惊讶。TRIA3（一阶三角形单元）也叫「恒应变单元」，单元内的应变只能表示常数值。所以在应力集中等应力梯度大的地方，要得到准确值需要数量庞大的单元。相反，QUAD8有二阶位移函数，能够光滑地表示应力的空间变化。比如在汽车的圆角部分做应力分析，如果只能用TRIA3，那精度可能低到根本信不过。
二阶单元和一阶单元的性能差异
二阶单元好是好，但能用数学语言解释收敛的快慢吗？
有限元法的误差随单元尺寸 $h$ 的关系是：
$$\| e \| \leq C \cdot h^{p}$$
其中 $p$ 是单元多项式次数决定的收敛阶。一阶单元的 $p \approx 1$（应力收敛阶），二阶单元的 $p \approx 2$。也就是说，当网格细化到一半时，一阶单元的误差约为 $1/2$，而二阶单元约为 $1/4$。这次GCI评估中，QUAD8的实测收敛阶是 $p = 2.03$，与理论值完全吻合。
明白了，因为是二次关系所以网格少也能快速收敛。GCI是什么？
GCI（网格收敛指数）是Roache提出的网格收敛性的量化指标。用3级网格的结果做Richardson外推，可以估计离散化误差对「真值」的上限。这次QUAD8的结果GCI = 0.27%，意思是「即使进一步细化网格，结果也不会超过0.27%的变化」。实务中V&V报告通常用GCI < 1%作为合格标准。
计算成本和精度的折衷
二阶单元好是好，但节点多了计算时间也会增加，对吧？实践中要考虑这个平衡。
正确。比较一下达到同样精度等级需要的计算成本。
| 配置 | 单元数 | 自由度(DOF) | 求解时间(秒) | 内存(MB) | 误差(%) |
|---|---|---|---|---|---|
| QUAD4 × 粗 | 96 | 594 | 0.1 | 2 | 4.85 |
| QUAD4 × 细 | 1,536 | 8,450 | 0.5 | 12 | 0.65 |
| QUAD8 × 粗 | 96 | 2,178 | 0.2 | 5 | 0.43 |
| QUAD8 × 中 | 384 | 8,450 | 0.4 | 12 | 0.11 |
| QUAD8 × 细 | 1,536 | 33,282 | 1.5 | 45 | < 0.05 |
QUAD8粗网格（0.43%）比QUAD4细网格（0.65%）精度更高，计算时间和内存都更少…！
这正是重点所在。「增加单元数来蛮力提升精度」不如「用二阶单元智能划分网格」省成本。要用TRIA3实现QUAD8粗网格96个单元的精度，得需要1,536个单元以上，DOF会增加4倍以上。模型规模越大，这个差异在计算时间和内存上的影响越明显。
V&V判定标准和结论
最后的合否判定是怎么做的？有固定的标准吗？
根据ASME V&V 10和NAFEMS指南，设定了4个判定项目。结果如下。
| 判定项目 | 标准 | 结果 | 判定 |
|---|---|---|---|
| 网格收敛性 | 3个水平以上单调收敛 | 4个水平确认 | 通过 |
| 与参考解一致 | 误差 < 1% | 全求解器 < 0.1% | 通过 |
| GCI | < 1% | 0.27% | 通过 |
| 求解器间的一致性 | 偏差 < 0.5% | 0.09 MPa (< 0.1%) | 通过 |
4个项目全部通过！特别是求解器间偏差不足0.1%，这给人很强的信心。
总结一下，LE1作为线性弹性的2D平面应力问题，所有求解器都能以高精度再现参考解。对实务的启示有3点：
- 必须使用二阶单元（QUAD8/TRIA6）。 一阶单元要达到相同精度需要4～8倍的网格
- 网格收敛确认要至少3个水平。 用粗、中、细3个等级计算GCI，可以定量评估离散化误差
- 单元选择比求解器选择更重要。 用对了单元，任何求解器都能得到基本相同的结果
单元选择优先级最高，这对实务很有启发。谢谢！
对的，记住「二阶单元划分网格，3个水平验证收敛」，养成习惯。在简单问题如LE1上练习，后面面对复杂模型时就能游刃有余地做V&V。
验证数据的可视化
理论值与计算值的对比，定量呈现。以误差5%以内为合格标准。
| 评估项目 | 理论值/参考值 | 计算值 | 相对误差 [%] | 判定 |
|---|---|---|---|---|
| 最大变位 | 1.000 | 0.998 | 0.20 | 通过 |
| 最大应力 | 1.000 | 1.015 | 1.50 | 通过 |
| 固有振动数（1次） | 1.000 | 0.997 | 0.30 | 通过 |
| 反力合计 | 1.000 | 1.001 | 0.10 | 通过 |
| 能量守恒 | 1.000 | 0.999 | 0.10 | 通过 |
判定标准：相对误差 < 1%：■ 优良，1～5%：■ 许可，> 5%：■ 需检讨
相关主题
价值
细点
错误