MMS：二维弹性体
理论：先定答案再检验代码
概述
怎样确认自己写的有限元代码是否正确？有精确解的问题很少吧。
这时就用制造解法（MMS）。先任意确定作为答案的位移，再反算使其成立的体力，这样任何形状都能构造出有精确解的问题。然后看网格加密时误差是否按理论速度减小。如果连速度（收敛阶）都一致，从方程组装、数值积分到边界条件，代码的大部分就可以说是正确的。我们用自编代码实际检验，并故意埋入缺陷看看会怎样。
制造解与体力
单位正方形，平面应力，弹性模量1，泊松比0.3。体力 $\mathbf{b}$ 用符号计算（sympy）推导，所有边界给定制造解位移。
理论收敛阶
$p$ 为单元阶数（Q4为1，Q8为2）。$p_{obs}$ 为网格加倍时由误差之比求得的观测阶数。
“制造答案”的逆向思维
制造解的思路从1980年代就在使用，但把它推广为代码验证方法的是Roache以及美国桑迪亚国家实验室的Salari和Knupp等人。他们在2000年前后的报告和论文中整理了用制造解检验收敛阶的步骤，指出“阶数一致就可以说不存在影响阶数的错误”。不是去找精确解而是去造精确解，这种逆向思维使方法可用于任何方程和代码。
算例：自编代码的收敛
正确的代码
| 单元·划分 | 位移误差 L2 | 阶数 | 梯度误差 H1 | 阶数 |
|---|---|---|---|---|
| Q4, 8 | 2.25e-3 | 2.04 | 7.74e-2 | 0.99 |
| Q4, 16 | 5.58e-4 | 2.01 | 3.88e-2 | 1.00 |
| Q4, 32 | 1.39e-4 | 2.00 | 1.94e-2 | 1.00 |
| Q4, 64 | 3.48e-5 | 2.00 | 9.70e-3 | 1.00 |
| Q8, 8 | 1.41e-4 | 2.93 | 7.37e-3 | 2.00 |
| Q8, 16 | 1.78e-5 | 2.99 | 1.85e-3 | 2.00 |
| Q8, 32 | 2.23e-6 | 3.00 | 4.62e-4 | 2.00 |
埋入缺陷的代码
| 单元·划分 | 位移误差 L2 | 阶数 |
|---|---|---|
| Q4·ν错误, 8 | 2.55e-3 | 1.88 |
| Q4·ν错误, 16 | 1.22e-3 | 1.06 |
| Q4·ν错误, 32 | 1.07e-3 | 0.19 |
| Q4·ν错误, 64 | 1.06e-3 | 0.02 |
| Q8·中点错误, 8 | 4.42e-4 | 2.23 |
| Q8·中点错误, 16 | 1.05e-4 | 2.08 |
| Q8·中点错误, 32 | 2.56e-5 | 2.03 |
泊松比的缺陷在8等分以前阶数为1.9，看起来很正常。
粗网格上离散误差更大，缺陷造成的小偏差被掩盖。网格加密、离散误差变小后，缺陷造成的1.1×10⁻³偏差原样留下，阶数降到0。所以至少要用四级网格确认阶数稳定在理论值。Q8中点的缺陷更难发现：误差一直很小，只是阶数从3降到2。因为边界值用直线连接，二次单元的能力在边界上没有发挥出来。这正是因为MMS看的是阶数而不是误差大小才能发现的缺陷。
MMS的做法
- 选择光滑、非平凡的制造解，使所有项（各应变分量、边界）都起作用。
- 体力用符号计算推导，避免手算错误。
- 网格按2倍加密四级以上，把观测阶数列表。
- 确认位移（L2）和梯度（H1）两者的阶数。
- 误差积分用足够点数的高斯积分，不混入积分误差。
“基准题都通过了”
某研究室的代码在几个基准题上与参考值相差1%以内，被认为是正确的。用MMS检查发现，应力误差的阶数比理论值低约0.5。原因是在弯曲边界上积分载荷时使用的点数不够。与基准参考值比较只能告诉你误差是否低于某个水平，而MMS的阶数能灵敏地显示某处存在错误。
常见错误
错误与对策
| 错误 | 后果 | 对策 |
|---|---|---|
| 只用两三级网格 | 缺陷被掩盖 | 四级以上 |
| 只看误差大小 | 漏掉阶数下降 | 计算观测阶数 |
| 制造解过于简单 | 部分项未被检验 | 所有分量都起作用的函数 |
| 手工推导体力 | 混入推导错误 | 用符号计算 |
| 误差积分点少 | 阶数算错 | 高阶高斯积分 |
我想了解相关内容。
相关页面参考MMS概述、MMS：二维导热、MMS收敛阶、代码验证和NAFEMS LE1。
相关课题
价值
更详细
错误