拉梅厚壁圆筒解 — 内压问题的FEM验证基准
拉梅厚壁圆筒解的理论基础
概述 — 为什么拉梅解最适合FEM验证
拉梅公式在压力容器教科书中出现过，但在FEM验证中也有用吗？
不仅有用，可以说是最优的轴对称基准。内压下厚壁圆筒的应力分布——半径方向应力 $\sigma_r$ 和周向应力 $\sigma_\theta$ ——的解析解已知且呈闭式，所以FEA数值解可以"一对一"进行答案检验。
与其他基准问题相比，比如悬臂梁弯曲，拉梅解有什么优势呢？
有三个优点：
- 应力梯度陡峭 — 应力在内表面达到最大值，向外表面急速衰减。这种条件下，网格质量和单元阶次的差异会明显体现出来
- 2次单元粗网格也精确 — 四边形8节点轴对称单元（QUAD8）只需4层网格，误差就可控制在0.12%以内。8层网格时误差实际为零
- 与工程实务直接相关 — ASME第VIII卷压力容器设计计算值与之直接对应。这意味着"教科书练习题"与"工程设计验证"可以用同一个模型完成
所以这不仅仅是理论问题，还能用于实际设计验证。NAFEMS基准集中有包括拉梅圆筒问题吗？
当然有。NAFEMS的LE10和LE11都属于厚壁圆筒类问题，ASME V&V 10-2006在推荐的Code Verification步骤中，也把"具有闭式解析解的问题"的代表列举出来，拉梅问题就在其中。部署新求解器或版本升级后回归测试时，这个模型应该是最先运行的问题之一。
控制方程 — 拉梅公式
请告诉我具体数学表达式。圆筒内径为 $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$ 和Hooke定律：
$\sigma_z$ 与半径位置 $r$ 无关，是常数。这也可以用来验证FEA结果，如果 $\sigma_z$ 在单元间波动，说明平面应变拘束设置有问题。
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 | 一般 |
| CAX4（4节点轴对称） | 8 | 100 | 165.5 | $-99.1$ | 0.72 | 通过 |
| CAX4（4节点轴对称） | 16 | 200 | 166.4 | $-99.8$ | 0.18 | 通过 |
| CAX8（8节点轴对称） | 4 | 100 | 166.5 | $-99.9$ | 0.12 | 通过 |
| CAX8（8节点轴对称） | 8 | 200 | 166.7 | $-100.0$ | 0.00 | 通过 |
| C3D20（20节点六面体3D） | 8×16×4 | 18,000 | 166.6 | $-99.9$ | 0.06 | 通过 |
判定标准：相对误差 < 1%: ■ 通过、1～5%: ■ 一般、> 5%: ■ 不通过
2次单元（CAX8）的收敛速度太快了！仅4分割就达到0.12%的误差，8分割时完全吻合。而1次单元就算4分割也要2.64%的误差。
这正是拉梅问题的教育价值。"为什么推荐2次单元"这个问题，从这个对比中就能直观理解。应力分布呈 $1/r^2$ 形式，有明显的曲率。1次单元的线性插值无法精确表示这种曲率，而2次单元的形状函数能直接表达这个曲率，所以用很少的单元就能达到高精度。
拉梅厚壁圆筒解的数值计算方法
轴对称模型的建立
用FEM解厚壁圆筒时，通常用轴对称单元，但也应该验证一下3D单元的结果吧？
从Code Verification角度，首先用轴对称2D求解。建模简单，自由度少，结果解释清晰。然后再用3D单元（例如1/4对称的90度扇形）验证求解器的3D模块能否给出相同结果——这样才能确保3D求解器正确处理了轴对称问题。
轴对称单元中，位移场只有两个分量 $u_r(r, z)$ 和 $u_z(r, z)$，应变向量是4分量：
其中 $\varepsilon_\theta = u_r/r$ 是轴对称特有的项，分母出现 $r$，所以在内表面附近数值上特别敏感。
单元刚性矩阵包含 $2\pi r$ 的体积分：
单元选择与积分方案
积分方法——完全积分和降阶积分——会改变结果吗？
在拉梅问题中影响不大，但了解差别很重要：
- 完全积分（2×2表示QUAD4，3×3表示QUAD8）：刚度偏大，但拉梅问题因为不是薄壳，不会出现严重的锁定现象
- 降阶积分（1×1表示QUAD4，2×2表示QUAD8）：计算效率更高，但1次单元存在沙漏模态风险。对拉梅问题，单元内位移模式简单，通常没有问题
- 非协调单元（Abaqus的CAX4I等）：在1次单元基础上添加额外位移模式以提高精度。如果一定要用1次单元做拉梅问题，推荐使用
| 单元 | 求解器名称 | 节点数 | 积分点 | 对拉梅问题的推荐度 |
|---|---|---|---|---|
| 4节点轴对称（完全积分） | CAX4 / PLANE182 / CQUAD4 | 4 | 2×2 | 基本验证 |
| 4节点轴对称（非协调模式） | CAX4I / PLANE182 w/ K-opt | 4 | 2×2 | 1次单元的首选 |
| 8节点轴对称（完全积分） | CAX8 / PLANE183 / CQUAD8 | 8 | 3×3 | 最推荐 |
| 8节点轴对称（降阶积分） | CAX8R / PLANE183 R | 8 | 2×2 | 推荐 |
网格收敛性与GCI
看了对比表能理解收敛趋势，但怎样定量判断"什么网格就足够了"呢？
用GCI（Grid Convergence Index，网格收敛指数）。这是Roache在1998年提出的方法，ASME V&V 20-2009也推荐使用。通过三个不同网格密度（粗、中、细）的结果，用Richardson外推法估计离散化误差的95%信赖区间。
具体地，对网格尺寸 $h_1 > h_2 > h_3$ 和对应的解 $f_1, f_2, f_3$：
其中 $r = h_1/h_2 = h_2/h_3$（等比网格细分）。观测收敛阶数 $p_\text{obs}$ 若接近理论值（1次单元：2，2次单元：3～4），说明网格处于渐近收敛区。
拉梅问题用2次单元时，通常半径方向4→8→16分割三水平会得到 $p_\text{obs} \approx 4$，GCI$_{12}$（8→16分割的信赖区间）小于0.01%。这说明8分割的网格已经充分，可以定量证明。
拉梅厚壁圆筒解的工程应用
分析流程与边界条件
在FEA软件中搭建拉梅问题时，从头到尾要做什么？
分步骤说明：
- 建立几何：在 $r$-$z$ 平面上创建矩形区域（$a \le r \le b$, $0 \le z \le L$）。$L$ 的值任意（因为平面应变，结果不受影响）
- 定义材料：线性弹性（$E = 210$ GPa, $\nu = 0.3$）
- 划分网格：半径方向用偏置网格（内表面密集）。推荐CAX8单元，至少4层
- 边界条件：
- $z = 0$ 面：$u_z = 0$（对称面）
- $z = L$ 面：$u_z = 0$（平面应变约束）
- $r = a$ 面：内压 $p = 100$ MPa（面荷载）
- $r = b$ 面：自由（外压 = 0）
- 求解：静态线性分析。用直接法求解器即可
- 后处理：沿半径方向提取 $\sigma_r$、$\sigma_\theta$、$\sigma_z$、$u_r$，与理论解叠图对比
边界条件中"$z = L$ 面设 $u_z = 0$"这条，如果漏掉会怎样？
漏掉的话就不再是平面应变，而变成广义平面应变状态，$\sigma_z$ 会改变，导致 $u_r$ 也偏离理论值5～15%。这是FEM验证中最常见的错误之一。
商用求解器的设置要点
Abaqus、Ansys、Nastran分别有什么要注意的地方吗？
三大主要求解器的设置对照：
| 项目 | Abaqus | Ansys Mechanical | MSC Nastran |
|---|---|---|---|
| 单元类型 | CAX8R（推荐）/ CAX4I | PLANE183（KEYOPT(3)=1） | CQUAD8 + PSHELL (AXI) |
| 内压荷载 | *DSLOAD, P1=100 | SFE, PRES=100 | PLOAD4 |
| 平面应变约束 | *BOUNDARY, ZSYMM | D, UZ, 0 on both ends | SPC1, UZ=0 |
| 应力输出 | S11=$\sigma_r$, S22=$\sigma_\theta$, S33=$\sigma_z$ | SX, SY, SZ（须确认坐标系） | 圆柱坐标系下的von Mises |
| 注意事项 | 从节点坐标提取$r$值 | 必须定义圆柱坐标系（LOCAL=11） | 用CORD2C定义圆柱坐标系 |
Ansys中如果坐标系设置错了，应力值会差得很远吧……
非常常见的错误。Ansys默认用直角坐标系（XYZ）输出应力，如果不转换到圆柱坐标系，$\sigma_r$ 和 $\sigma_\theta$ 会完全错误。Abaqus的轴对称单元从一开始就是 S11=$r$ 方向、S22=$\theta$ 方向，所以较少出现这种混淆。但3D模型中同样会有坐标系问题。
ASME第VIII卷的一致性检验
拉梅公式在ASME压力容器规范中也用到吧。设计值与FEA对比时有什么要点吗？
ASME第VIII卷第2版规定，厚壁圆筒的许用内压为：
这里 $S$ 是许用应力，$k = b/a$。这个式子实际上是从拉梅公式反推得出的，令 $\sigma_\theta(a) = S$ 就得到这个式子。
工程中要注意的点：
- Division 1 vs Division 2：Division 1采用薄壳近似简式，Division 2采用拉梅严格解。FEA验证应该对标Division 2
- 腐蚀余量（CA）：ASME计算时要减去腐蚀余量的有效厚度。FEA验证用的是公称尺寸（无余量）
- 闭端效应：实际压力容器有鏡板，轴向产生 $\sigma_z = pa^2/(b^2 - a^2)$ 的恒定应力，与平面应变的 $\sigma_z$ 值不同。需要对齐验证条件
拉梅厚壁圆筒解的软件对比
求解器结果对比
用多个求解器解同一个问题，结果会有差异吗？
对拉梅问题这样的线性弹性轴对称问题，用相同的单元和网格，求解器间的差异基本为零。反过来说，如果求解器间有差异，那就说明某处设置有误，这正是拉梅问题的价值所在。
| 求解器 | 单元 | 半径方向8分割 | $\sigma_\theta(a)$ [MPa] | 误差 [%] |
|---|---|---|---|---|
| Abaqus 2024 | CAX8R | 8 | 166.67 | 0.00 |
| Ansys 2024R2 | PLANE183 | 8 | 166.67 | 0.00 |
| MSC Nastran 2023 | CQUAD8 (AXI) | 8 | 166.66 | 0.01 |
| COMSOL 6.2 | 二阶三角形 | 8相当 | 166.65 | 0.01 |
| CalculiX 2.21 | CAX8 | 8 | 166.67 | 0.00 |
| Code_Aster 15.8 | QUAD8_AXI | 8 | 166.66 | 0.01 |
开源软件CalculiX和Code_Aster的精度也和商用软件一样好。这样的话，新求解器评估时能不能用拉梅问题验证？
完全可以。"先用拉梅问题验证新求解器，看6位数是否一致"是业界的标准做法。倒过来说，拉梅问题中误差很大的求解器，其基础实现可能有问题，应当谨慎使用。
拉梅厚壁圆筒解的前沿研究
扩展问题 — 热应力·弹塑性·蠕变
拉梅的线性弹性解理解了。但实际的压力容器工作在高温，还要降伏，有时还出现蠕变。这些情况也能用拉梅问题作为基准吗？
问得很好。拉梅问题的推广形式在以下方向都有解析解或半解析解，可以作为更高级的验证基准：
- 热应力（NAFEMS LE11）：内外表面加温度梯度，内表面应力叠加热应力。温度呈对数分布，应力呈 $1/r^2$ 分布，网格对温度和应力分布的同时捕捉能力要求高，比单纯的力学问题更能暴露网格细节问题
- 弹塑性（Hill解）：逐步增加内压，塑性区从内表面向外扩展。弹塑性边界半径 $c$ 有理论公式。这个问题同时考查线性求解和非线性算法
- 稳态蠕变（Norton律）：在蠕变应力律 $\dot{\varepsilon} = A\sigma^n$ 下，应力重新分布，内表面的集中逐渐平缓化。这样可以验证求解器的蠕变本构模型
- 自紧处理（Autofrettage）：预先加高压使内表面塑性变形，然后卸载残留压应力。结合疲劳寿命计算，可定量评估自紧的效果
同一个几何，从线性弹性一路验证到弹塑性、热连成、蠕变……这样逐步增加复杂性的方法很有效啊。
正是这样。这叫"验证阶梯"（Verification Ladder）。对FEM初学者，先在线性弹性拉梅问题上完全掌握，再逐步扩展到弹塑性和热连成，学习曲线会平缓很多。每一步都有理论解作为参照，能够明确看到"新增模块改变了什么"。
拉梅厚壁圆筒解的故障排除
FEM验证的常见错误
"理论解和FEA不吻合！"这种情况下，通常是什么原因导致的？
根据过往经验，常见失败原因排名如下：
| 排名 | 失败模式 | 表现 | 对策 |
|---|---|---|---|
| 1 | 平面应变约束缺失 | $\sigma_z$ 变成零，$u_r$ 离理论值5～15% | 检查z方向两端是否都设了UZ=0。有些求解器有"广义平面应变"的设置选项 |
| 2 | 坐标系混淆 | $\sigma_r$ 和 $\sigma_\theta$ 对调，或用直角坐标成分对比 | 确认轴对称单元应力分量定义（Abaqus: S11=$r$, S22=$\theta$）。3D模型则必须定义圆柱坐标系并变换 |
| 3 | 荷载方向或面选择错误 | 总体符号相反，或外表面受压而非内表面 | 查看荷载方向定义（内向vs外向）和位置。用变形图目视检查膨胀方向是否正确 |
| 4 | 单位制不一致 | 应力数值大小差很多倍 | 统一成 mm-N-MPa 系或 m-N-Pa 系。检查 $E = 210,000$ MPa = $2.1 \times 10^{11}$ Pa |
| 5 | 1次单元精度极限 | "细分网格到32层也降不到2%误差" | CAX4细分不如CAX8粗分。建议先用CAX8搭4分割确认"正确答案"，再研究CAX4的限制 |
| 6 | 应力外推值vs积分点值 | 节点外推值和积分点值应力不同（特别是内表面） | 拉梅问题验证用积分点值（单元输出）。节点外推会因平均化而钝化，特别在应力梯度大的区域 |
排名第1的"平面应变约束缺失"这种基本问题也容易出现啊。自己一定会犯。
越是基本的东西越容易被忽视。应对办法是分析完就立刻检查 $\sigma_z$ 的值。理论值是 $\sigma_z = 2\nu a^2 p/(b^2 - a^2)$。本例中 $\sigma_z = 20.0$ MPa，而且要在所有单元中都是这个常数。如果 $\sigma_z \approx 0$，说明变成平面应力了；如果不均匀波动，说明约束有问题——立刻就能判断。
$\sigma_z$ 当作"金丝雀"来用，我记住了！有没有更完整的排查清单？
最后总结一套调试流程，遇到"理论解和FEA不吻合"时照着来：
- 看变形图 — 膨胀方向正确吗？有异常的变形模式吗？
- 检查 $\sigma_r + \sigma_\theta$ 是否为常数 — 拉梅解的基本性质。波动就说明设置有问题
- 检查 $\sigma_z$ 是否为常数 — 平面应变约束的验证
- 检查外表面 $\sigma_r = 0$ — 外压零的边界条件
- 检查内表面 $\sigma_r = -p$ — 内压的大小和符号
这5项全部通过，就基本能确保与理论值一致。
详情
报告