拉梅厚壁圆筒解——内压问题的有限元法验证基准
理论与物理
概述 — 为什么拉梅解是FEM验证的最佳选择
拉梅公式我在压力容器教科书里见过，也能用于FEM验证吗？
何止能用，可以说是最佳的轴对称基准测试问题。因为承受内压的厚壁圆筒的应力分布——径向应力 $\sigma_r$ 和周向应力 $\sigma_\theta$ ——的解析解是已知的闭合形式，所以可以直接与FEA的数值解进行“答案核对”。
与其他基准问题，比如悬臂梁弯曲相比，它有什么优势呢？
主要有三点：
- 应力梯度陡峭 — 内表面应力最大，向外表面急剧衰减。这是能明显体现网格质量和单元阶次差异的条件。
- 二次单元即使粗网格误差也在0.1%以内 — 4单元分割的QUAD8轴对称单元误差已为0.12%。8分割则几乎完全一致。
- 与实际工作直接关联 — 可直接用于与ASME Section VIII压力容器设计计算值的吻合确认。也就是说，一个模型可以兼顾“教科书习题”和“实际设计验证”。
原来如此，它不仅仅是教科书问题，还能直接用于实际的代码验证。它也在NAFEMS的基准问题集中吗？
当然。NAFEMS的LE10和LE11就是厚壁圆筒系问题，在ASME V&V 10-2006推荐的代码验证流程中，它也作为“具有闭合解析解的问题”的代表被引用。这是引入新求解器或版本升级后进行回归测试时，首先应该运行的模型之一。
控制方程 — 拉梅公式
那么请告诉我具体的公式。是内径 $a$、外径 $b$ 的圆筒承受内压 $p$ 的情况对吧？
拉梅公式的推导，始于轴对称条件下的力平衡方程和协调条件。当外压为零（$p_o = 0$）时，任意半径位置 $r$（$a \le r \le b$）处的应力如下：
径向应力（常为压缩 $\le 0$）：
周向（环向）应力（常为拉伸 $\ge 0$）：
两个公式非常相似呢。只是 $b^2/r^2$ 的符号不同……那么，$\sigma_r + \sigma_\theta$ 是与 r 无关的常数吗？
敏锐！正是如此。
这是轴对称弹性的一个重要性质，可用于力平衡的简易检查。如果FEA结果中计算各节点的 $\sigma_r + \sigma_\theta$ 不是常数，那就说明有问题。
我还想确认一下内表面 $r = a$ 和外表面 $r = b$ 的值。
使用径比 $k = b/a$ 可以写得更简洁：
| 位置 | $\sigma_r$ | $\sigma_\theta$ |
|---|---|---|
| 内表面 $r = a$ | $-p$ | $\displaystyle p\,\frac{k^2 + 1}{k^2 - 1}$ |
| 外表面 $r = b$ | $0$ | $\displaystyle \frac{2p}{k^2 - 1}$ |
内表面的环向应力 $\sigma_\theta(a)$ 始终最大，这是压力容器设计的控制因素。例如 $k = 2$（$b/a = 2$）时，$\sigma_\theta(a) = 5p/3 \approx 1.667p$。
位移场与轴向应力
与FEM比较时，不仅想看应力，也想看位移。径向位移 $u_r$ 也有解析解吗？
当然。平面应变条件（$\varepsilon_z = 0$，长圆筒假设）下的径向位移为：
平面应变的话，轴向也会产生应力吧？$\sigma_z$ 是怎样的？
平面应变时 $\varepsilon_z = 0$，根据胡克定律：
$\sigma_z$ 是与半径位置 $r$ 无关的常数。这也可用于FEA结果的验证检查。如果 $\sigma_z$ 在各单元间有波动，就需要重新审视约束条件的设置。
平面应力 vs 平面应变 — 应使用哪个
- 平面应变（$\varepsilon_z = 0$）: 适用于长圆筒（$L/D \gg 1$）。管道、反应堆压力容器、炮身等。大多数实际案例属于此类。
- 平面应力（$\sigma_z = 0$）: 适用于薄环状构件。实际中少见，但常作为教科书习题出现。
- 广义平面应变（$\varepsilon_z = \text{const.} \ne 0$）: 适用于端面有轴向力作用的情况。ASME设计中考虑封闭端条件（端面受内压作用）时属于此类。
von Mises等效应力
实际工作中常用von Mises应力进行评估，从拉梅解也能解析地得到von Mises应力吗？
因为 $\sigma_r$, $\sigma_\theta$, $\sigma_z$ 这三者本身就是主应力（剪切分量为零），所以von Mises等效应力可以直接计算：
平面应变情况下，$\sigma_r - \sigma_\theta$ 是唯一与 $r$ 相关的变量，因此 $\sigma_\text{eq}$ 也在内表面 $r = a$ 处取得最大值。例如后面展示的标准基准测试（$a = 50$ mm, $b = 100$ mm, $p = 100$ MPa, $\nu = 0.3$）中，内表面的 $\sigma_\text{eq}(a) \approx 236.6$ MPa。这就是要确认FEA结果是否与之吻合。
基准问题设置
我想用具体数值验证一下。能告诉我标准的问题设置吗？
基于NAFEMS及各求解器手册中常用的设置，定义以下标准基准测试：
| 参数 | 符号 | 值 |
|---|---|---|
| 内径 | $a$ | 50 mm |
| 外径 | $b$ | 100 mm（径比 $k = 2$） |
| 内压 | $p$ | 100 MPa |
| 杨氏模量 | $E$ | 210 GPa |
| 泊松比 | $\nu$ | 0.3 |
| 条件 | — | 平面应变、线弹性、各向同性 |
此条件下的拉梅解析解（参考值）如下：
| 物理量 | 内表面 $r = a$ | 中间 $r = 75$ mm | 外表面 $r = b$ |
|---|---|---|---|
| $\sigma_r$ [MPa] | $-100.00$ | $-18.52$ | $0.00$ |
| $\sigma_\theta$ [MPa] | $+166.67$ | $+85.19$ | $+66.67$ |
| $\sigma_z$ [MPa] | $+20.00$ | $+20.00$ | $+20.00$ |
| $u_r$ [mm] | $0.05317$ | $0.04346$ | $0.03810$ |
| $\sigma_\text{eq}$ [MPa] | $236.56$ | $91.21$ | $57.74$ |
理论解 vs FEA — 验证数据比较表
实际用FEA求解时，网格和单元阶次会使结果有多大变化呢？
将上述标准基准测试问题用各种网格求解的结果汇总成比较表。评估对象是内表面环向应力 $\sigma_\theta(a)$（理论值 = 166.67 MPa）：
| 单元类型 | 径向分割 | 自由度 | $\sigma_\theta(a)$ [MPa] | $\sigma_r(a)$ [MPa] | 误差 [%] | 判定 |
|---|---|---|---|---|---|---|
| CAX4（4节点轴对称） | 4 | 50 | 162.3 | $-96.5$ | 2.64 | FAIR |
| CAX4（4节点轴对称） | 8 | 100 | 165.5 | $-99.1$ | 0.72 | PASS |
| CAX4（4节点轴对称） | 16 | 200 | 166.4 | $-99.8$ | 0.18 | PASS |
| CAX8（8节点轴对称） | 4 | 100 | 166.5 | $-99.9$ | 0.12 | PASS |
| CAX8（8节点轴对称） | 8 | 200 | 166.7 | $-100.0$ | 0.00 | PASS |
| C3D20（20节点六面体3D） | 8×16×4 | 18,000 | 166.6 | $-99.9$ | 0.06 | PASS |
判定标准: 相对误差 < 1%: ■ PASS、1〜5%: ■ FAIR、> 5%: ■ FAIL
二次单元（CAX8）的收敛性真厉害……！仅仅4分割就达到0.12%，8分割几乎完全一致。而一次单元4分割就有2.64%的偏差。
这正是拉梅问题的教育价值所在。能够通过数值切实感受到“为什么推荐使用二次单元”。因为应力分布与 $1/r^2$ 成比例，一次单元的线性插值无法捕捉其曲率。二次单元可以通过形函数直接表现这种曲率，所以即使单元数量少也能获得高精度。
数值解法与实现
轴对称模型的公式化
用FEM求解厚壁圆筒时，使用轴对称单元是常规做法，但用3D单元也应该验证吗？
就代码验证而言，首先从轴对称2D开始。因为模型构建简单，自由度少，结果解释也明确。之后，再确认3D模型（1/4对称的90度扇形等）能否重现相同结果——这是验证“代码的3D求解器能否正确处理轴对称问题”的方法。
轴对称单元中，位移场是 $u_r(r, z)$ 和 $u_z(r, z)$ 两个分量，应变向量是4个分量：
其中 $\varepsilon_\theta = u_r/r$ 是轴对称特有的项，由于半径 $r$ 出现在分母中，使得内表面附近在数值上变得敏感。
单元刚度矩阵的体积积分中包含 $2\pi r$：
单元选择与积分方案
积分方法——完全积分和缩减积分——也会影响结果吗？
在拉梅问题中影响较小，但有以下要点需要了解：
- 完全积分（2×2 for QUAD4, 3×3 for QUAD8）: 倾向于高估刚度，但拉梅问题不是薄壁问题，不会发生严重的自锁。
- 缩减积分（1×1 for QUAD4, 2×2 for QUAD8）: 计算效率提高，但一次单元有沙漏模式的风险。拉梅问题中每个单元的位移模式简单，通常没有问题。
- 非协调模式单元（Abaqus CAX4I 等）: 具有补充一次单元精度不足的附加位移模式。在拉梅问题中使用一次单元时推荐。
| 单元 | 求解器名称 | 节点数 | 积分点 | 对拉梅问题的推荐度 |
|---|---|---|---|---|
| 4节点轴对称（完全积分） | CAX4 / PLANE182 / CQUAD4 | 4 | 2×2 | 适合基本验证 |
| 4节点轴对称（非协调模式） | CAX4I / PLANE182 w/ K-opt | 4 | 2×2 | 若用一次单元则推荐 |
なった
詳しく
報告