矩形板的Navier解（四边简支）
理论：双重傅里叶级数精确解
概述
我建了一个壳单元的FEM模型，有办法确认结果是否正确吗？
对于板的弯曲，先与Navier解比较是标准做法。四边简支矩形板受均布荷载的问题，纳维于1820年用双重傅里叶级数给出了精确解。答案明确，可以一次性检验单元类型、网格疏密和边界条件是否正确，非常适合在进入实际模型之前“试手”。
控制方程与解
这是基尔霍夫薄板理论，$D$ 为弯曲刚度，$h$ 为板厚。边长 $a \times b$、四边简支、受均布荷载 $q_0$ 时，正弦函数之积自动满足边界条件，解为如下级数（$m, n$ 仅取奇数）：
系数表（泊松比0.3）
| 边长比 b/a | 挠度系数 α | 弯矩系数 β（短边方向） |
|---|---|---|
| 1.0 | 0.00406 | 0.0479 |
| 1.5 | 0.00772 | 0.0812 |
| 2.0 | 0.01013 | 0.1017 |
| 3.0 | 0.01223 | 0.1189 |
$a$ 为短边。边长比增大时，板趋近于沿短边方向跨越的梁（$\alpha = 5/384 \approx 0.01302$，$\beta = 1/8$）。
纳维曾想建造一座桥
克劳德-路易·纳维以流体的纳维–斯托克斯方程闻名，但他的本行是桥梁工程师，1820年对板的分析正是这一工作的延伸。然而他在巴黎塞纳河上规划的悬索桥，在1826年接近完工时因基础问题而被拆除。这常被用来说明理论的精确与实物的成功是两回事。解析解是检验“模型是否正确建立”的工具，并不保证“设计是否安全”。
算例：级数收敛与网格收敛
取多少项才够
| 所取项（m, n 上限） | 挠度系数 α | 弯矩系数 β |
|---|---|---|
| 1 | 0.004161 (+2.4 %) | 0.05338 (+11.5 %) |
| 3 | 0.004055 | 0.04692 |
| 5 | 0.004064 | 0.04823 |
| 9 | 0.004062 | 0.04797 |
| 19 | 0.004062 | 0.04788 |
| 收敛值 | 0.004062 | 0.04789 |
挠度大约取到5就定了，弯矩却迟迟不稳定。
弯矩是挠度的二阶导数，每一项都乘上 $m^2$ 或 $n^2$，高阶项的贡献更大。剪力是三阶导数，收敛更慢。FEM 也是如此：同一套网格，挠度已经吻合而应力还没收敛的情况很常见。验证时一定要同时比较“挠度”和“应力”。
钢板算例
1 m × 1 m、厚10 mm 的钢板（$E = 206$ GPa，$\nu = 0.3$）受5 kPa 均布荷载时，$D \approx 18{,}900$ N·m，中心挠度约1.08 mm（板厚的11%），最大弯曲应力约14 MPa。挠度相对板厚足够小，处于线性薄板理论的适用范围内。
差分法网格收敛示例
简支板可以拆分为两个泊松方程求解。用 $N \times N$ 网格的差分法求解正方形板，并把中心挠度系数与收敛值比较：
| 网格 | 挠度系数 α | 误差 |
|---|---|---|
| 4 × 4 | 0.004028 | −0.84 % |
| 8 × 8 | 0.004055 | −0.19 % |
| 16 × 16 | 0.004061 | −0.04 % |
| 32 × 32 | 0.004062 | −0.01 % |
网格每加密一倍误差约减为1/4，符合二阶方法的预期收敛。FEM 也一样，用三套以上网格确认误差的减小方式，就能判断单元是否正常工作。
与FEM比较时的注意事项
简支约束的施加
- 沿边缘只约束面外位移 $w = 0$。若连边缘转动也约束，就变成固支，挠度会减小到约0.3倍。
- 面内位移只做防止刚体移动的最小约束。若约束整条边的面内位移，挠曲时会产生膜力，尤其在大挠度时结果会偏离。
- Navier 解中四角会产生向下的集中反力，若无支承，角部会翘起。若FEM把支承建模为只受压的接触，角部会翘起，结果就与解析解不一致。
单元与板厚的选择
- 剪切锁定：完全积分的一次单元在薄板中过于刚硬，严重低估挠度。应使用减缩积分单元或具有防锁定措施的单元。
- 厚板：Navier解基于薄板理论，板厚超过边长约1/10时，FEM（Mindlin理论）会因剪切变形而得到更大的挠度。比较时应使用薄板模型。
- 大挠度：挠度达到板厚的相当比例后，膜效应会使实际挠度小于线性解。应线性对线性比较，非线性效应另行评估。
验证步骤
- 建立板厚约为边长1/100的正方形板，四边只约束 $w = 0$，施加均布荷载。
- 用三套以上网格，把中心挠度和弯矩（或表面应力）与系数表比较。
- 确认挠度和应力都收敛到1～2%以内后，再把同样的单元类型和网格密度用于实际模型。
“应力对不上”多半是读取位置的问题
壳单元的应力通常在单元中心或积分点计算，节点值多由周围单元平均得到。若板中心恰好没有节点，或直接比较节点平均后的云图，同一模型也会相差几个百分点。系数表中的 $\beta$ 是板中心、板表面的值，比较时要统一位置（中心）、表面（上表面或下表面）以及显示的平均方式。
常见错误
错误与对策
| 错误 | 后果 | 对策 |
|---|---|---|
| 同时约束边缘转动 | 变成固支，挠度约为0.3倍 | 只约束 $w = 0$ |
| 完全积分一次单元 | 剪切锁定导致过刚 | 减缩积分或防锁定单元 |
| 只凭挠度判断收敛 | 应力仍未收敛 | 应力也用三套网格确认 |
| 用厚板比较 | 因剪切变形而偏离 | 板厚取边长的1/100左右 |
| 使用不同泊松比的系数表 | 弯矩相差几个百分点 | 按材料的 $\nu$ 计算系数 |
还想了解其他验证用的问题。
轴对称问题可以用圆板弯曲，一维问题可以用简支梁等分布荷载。板的屈曲可用薄板屈曲计算器，层合板刚度可用复合材料层合板模拟器试算。