代码验证(Code Verification)
理论:用数字证明程序正确
概述
被要求证明公司自研分析程序的质量。用例题说“结果看起来合理”还不够吗?
不够。小错误只让结果稍有偏离,肉眼几乎看不出来。代码验证是在答案已知的问题上,确认误差能否小到计算机舍入误差量级,或网格加密时能否按理论阶数减小。有限元的基本方法是分片试验:单元应能表示的解,在单元畸变时能否精确再现。关键是合格与否能用数字清楚地表示出来。
试验用问题
两者都是拉普拉斯方程的精确解。一次式可被一次三角形单元精确表示,用于分片试验;指数与正弦之积用于阶数试验。
单元刚度与误差度量
$A^e$ 为三角形面积,$B$ 为形函数梯度矩阵。看节点最大误差以及网格间距减半时的阶数 $p$。
分片试验的诞生
分片试验是20世纪60年代中期Irons等人在开发板弯曲单元时提出的。当时有些单元计算能正常运行,但网格加密后却不趋近正确答案。Irons提出,对由几个单元组成的小“片”施加单元应能表示的简单状态(如均匀应变),以能否精确再现来判断单元是否正确。后来Strang等人阐明了其数学意义,至今仍是检验新单元和自研程序的第一个试验。
算例:小型有限元程序
例1:分片试验
单位正方形划分为6×6的一次三角形单元,内部节点最多偏移网格间距的30%使其畸变。边界给定 1 + 2x + 3y 求解,节点最大误差为1.8×10⁻¹⁵(舍入误差量级),证明实现正确。
例2:阶数试验与错误检测
用指数与正弦之积的解加密网格,比较正确实现与植入“x > 0.5范围单元面积算成一半”错误的实现:
| 分段数 | 正确实现误差 | 阶数 | 错误实现误差 | 阶数 |
|---|---|---|---|---|
| 4 | 5.8×10⁻⁴ | — | 8.8×10⁻² | — |
| 8 | 1.6×10⁻⁴ | 1.86 | 1.0×10⁻¹ | −0.21 |
| 16 | 4.1×10⁻⁵ | 1.98 | 1.0×10⁻¹ | 0.01 |
| 32 | 1.0×10⁻⁵ | 2.00 | 9.9×10⁻² | 0.02 |
错误实现的误差加密网格也完全不减小。
对。面积减半使该范围单元刚度翻倍,成为加密网格也修正不了的“方程本身的偏差”,所以阶数为0。同一个错误也使分片试验误差跳到0.23。反过来说,通过分片试验和阶数试验,就能用数字说明至少不存在这类错误。重要的是把这些试验自动化,每次修改程序都必须运行。一次通过不代表以后的修改不会把它弄坏。
代码验证步骤
- 用分片试验(均匀应变、线性温度分布等)确认单元基本性能。
- 在有精确解的问题上,确认加密网格能得到理论阶数。
- 没有精确解的方程或项,用构造解方法确认。
- 按边界条件类型、单元类型、非线性设置分别准备试验。
- 把试验自动化,每次修改程序都运行并记录。
“新版本结果变了0.3%”
某公司自研结构分析程序改进后的新版本,用同一模型计算,应力只差0.3%。差点被当作误差范围放过,但自动化的分片试验有一项失败了。调查发现,为提速而改写的单元积分只在畸变单元时把权重弄错了。0.3%看起来很小,但通不过分片试验的错误加密网格也不会消失。多亏有用数字判定合格与否的试验,才在发布前发现了问题。
常见错误
错误与对策
| 错误 | 后果 | 对策 |
|---|---|---|
| 凭看似合理的结果判断 | 漏掉小错误 | 用数字判定合否 |
| 只用规整网格试验 | 残留畸变单元错误 | 畸变的片 |
| 只用一套网格比较 | 漏掉阶数下降 | 加密求阶数 |
| 不试验边界条件 | 边界实现错误 | 按类型试验 |
| 只试验一次就结束 | 以后的修改使其损坏 | 自动化并每次运行 |
我想了解相关内容。
可以用Richardson外推工具试算阶数。相关页面参考构造解方法、解验证、Roache的GCI、ASME V&V 10和NAFEMS LE1。