简支梁等分布荷载 — FEM验证基准测试
简支梁等分布荷载的理论基础
概要
什么是简支梁基准测试？要检查什么？
这是FEM最基础的验证问题。承受等分布荷载 $w$ 的简支梁拥有来自Euler-Bernoulli梁理论的严格解析解。中央的最大挠度为
通过与数值解的比较，可以检验要素类型的精度。特别重要的是，采用一次要素（线性要素）时会出现剪切自锁现象，导致梁的表现过于刚硬，挠度被严重低估。定量把握这一现象，以及通过低阶积分或B-bar法改善的效果，是此基准测试的实务关键。
只有一个挠度公式就能说明这么多问题吗？
是的。正因为它够简单，所以用处才大。NAFEMS的基准测试集合中有收录，ASME V&V 10-2006的代码验证（Code Verification）的第一步也推荐这个问题。在引入新求解器时，首先要确认这个问题的结果与理论解在0.1%以内，这是行业的标准做法。比如汽车制造商的结构分析部门，每次求解器版本升级时都要把这个基准测试作为回归测试跑一遍。
支配方程
请告诉我完整的理论式。不仅是挠度，应力也要出。
从Euler-Bernoulli梁的支配方程开始。承受等分布荷载 $w$（N/m）、长度为 $L$ 的简支梁的挠度曲线为
其中 $x$ 是从左端起的距离，$E$ 是杨氏模量，$I$ 是截面惯性矩。中央 $x = L/2$ 代入得最大挠度。
接下来是弯矩分布。
中央处最大，$M_{\max} = wL^2/8$。从而最大弯曲应力为
$Z = I/c$ 为截面系数，$c$ 为中性轴到最外缘的距离。对于矩形断面 $b \times h$，有 $I = bh^3/12$，$c = h/2$，$Z = bh^2/6$。
还有剪力分布吧？
剪力为
支点处最大 $V_{\max} = wL/2$，中央处为零。Euler-Bernoulli理论基于截面保持平面且法线保持垂直的假设（忽略剪切变形），对于跨径/高度比 $L/h$ 较大的薄梁适用性高。
Timoshenko梁与比较
当 $L/h$ 很小时会怎样？与Timoshenko理论有什么区别？
Timoshenko梁理论考虑了剪切变形。中央挠度为
第二项是剪切变形导致的额外挠度。$\kappa$ 是剪切修正系数（矩形断面为 $5/6$），$G$ 是剪切弹性模量 $G = E/\{2(1+\nu)\}$。
| $L/h$ | 弯曲挠度比率 | 剪切挠度比率 | Euler-Bernoulli误差 |
|---|---|---|---|
| 50 | 99.7% | 0.3% | 0.3% |
| 20 | 98.1% | 1.9% | 1.9% |
| 10 | 92.9% | 7.1% | 7.1% |
| 5 | 78.1% | 21.9% | 21.9% |
当 $L/h \geq 20$ 时Euler-Bernoulli理论足够准确。$L/h < 10$ 的深梁需要Timoshenko理论或3D固体要素。基准测试通常采用 $L/h = 50$ 左右，以确保Euler-Bernoulli理论严格成立的条件下进行验证。
基准验证数据
可以给个具体的数值例子吗？标准的参数设置是什么？
标准基准测试设置如下。
| 参数 | 记号 | 数值 |
|---|---|---|
| 跨径 | $L$ | 1000 mm |
| 截面宽 | $b$ | 100 mm |
| 截面高 | $h$ | 20 mm |
| 等分布荷载 | $w$ | 1.0 N/mm |
| 杨氏模量 | $E$ | 200,000 MPa |
| 泊松比 | $\nu$ | 0.3 |
$L/h = 50$ 满足Euler-Bernoulli理论条件。截面惯性矩 $I = bh^3/12 = 100 \times 20^3/12 = 66{,}667$ mm$^4$。
其中 $Z = bh^2/6 = 100 \times 20^2/6 = 6{,}667$ mm$^3$。
那不同要素类型的结果如何？
用上述参数的比对结果如下。
| 要素类型 | 网格 | 自由度 | $\delta_{\max}$ [mm] | 挠度误差 [%] | $\sigma_{\max}$ [MPa] | 应力误差 [%] |
|---|---|---|---|---|---|---|
| BEAM2（Euler-Bernoulli） | 20单元 | 126 | 0.9766 | 0.00 | 18.75 | 0.00 |
| BEAM3（Timoshenko） | 20单元 | 126 | 0.9795 | 0.30 | 18.75 | 0.00 |
| QUAD8（二次壳） | 20×2 | 2,520 | 0.9764 | 0.02 | 18.72 | 0.16 |
| HEX8 完全积分 | 50×10×4 | 19,800 | 0.8112 | 16.93 | 17.92 | 4.43 |
| HEX8 低阶积分 | 50×10×4 | 19,800 | 0.9738 | 0.29 | 18.70 | 0.27 |
| HEX20（二次固体） | 20×4×2 | 15,120 | 0.9762 | 0.04 | 18.73 | 0.11 |
| TET4（线性四面体） | 自动 | ~30,000 | 0.6834 | 30.02 | 15.21 | 18.88 |
| TET10（二次四面体） | 自动 | ~25,000 | 0.9741 | 0.26 | 18.68 | 0.37 |
HEX8完全积分产生约17%的挠度误差，TET4更是高达30%！这就是剪切自锁的威力。而低阶积分的HEX8和HEX20则误差都在0.3%以内。简支梁基准测试的意义就在于清晰地展示这一差异。
简支梁等分布荷载的数值计算方法
剪切自锁的物理
刚才表格中HEX8和TET4的结果这么差，是「剪切自锁」的缘故吧？具体发生了什么？
要理解本质，需要看一次要素形函数的局限性。
一次六面体要素（HEX8）的位移场在各方向都是线性的。纯弯曲时，截面会转动，上下表面产生相反方向的水平位移。理论上这时剪切应变 $\gamma_{xy} = 0$ 才对，但线性位移场无法满足这个条件。
也就是说，出现了本不应存在的寄生剪切应变（parasitic shear）。这个虚假的剪切应变对应着剪切应力，产生额外的刚性，使得梁的实际挠度被大幅低估。这就是剪切自锁。
明白了……要素的形函数没办法表达弯曲变形模式，所以才会产生「假剪切」。
正是。细化网格能缓解，但收敛太慢。厚度方向只有2个单元远远不够，需要8～16个单元才行，计算成本会暴增。所以实务上，采用自锁对策技术要高效得多。
自锁对策技术
有什么办法可以对付剪切自锁？
主要有三个办法。
1. 低阶积分（Reduced Integration）
不用完全积分（HEX8为2×2×2=8点），改用低阶积分（1×1×1=1点）。积分点数减少，寄生剪切应变的能量贡献被排除了。但代价是出现沙漏模式（hourglass mode）——零能量的变形模式，看起来像砂漏被捏变形，但计算中能量为零。
2. B-bar法（选择性降阶积分）
体积变化部分用低阶积分，偏差部分用完全积分。既能消除体积自锁又能消除剪切自锁，两个问题都解决。Abaqus的C3D8R配合沙漏控制，Nastran的CHEXA降阶版本，都是这个思路。
3. EAS法（增强假设应变法）
在标准位移场基础上，要素内部增加额外的应变模式。内部自由度通过静止凝聚对外界不可见。Ansys的SOLID185（keyopt(2)=2）用的就是这个。
| 方法 | 优点 | 注意事项 |
|---|---|---|
| 低阶积分 | 实现简单，计算成本低 | 必须配合沙漏控制 |
| B-bar法 | 体积和剪切自锁双解决 | 需要对要素技术有深入理解 |
| EAS法 | 精度高，无沙漏问题 | 计算成本稍高，非线性时偶有不稳定 |
| 二次要素 | 根本解决，精确表达弯曲模式 | DOF数约3倍，网格生成较费力 |
沙漏模式的控制具体怎么做？
沙漏模式的抑制需要加上人为的刚性（hourglass stiffness）。如果沙漏能量占总内部能的5～10%以上，就是危险信号。Abaqus用 *SECTION CONTROLS, HOURGLASS=ENHANCED，LS-DYNA用 IHQ=6（Belytschko-Bindeman）都是经验丰富的设置。这个基准问题中，如果沙漏能量占比在1%以下，就可以判定为合格。
网格收敛性理论
网格细化时，结果离理论解的距离下降有多快？
根据FEM的先验误差估计（Cea定理），$p$ 阶要素的能量范数误差以 $O(h^p)$ 下降，$L^2$ 范数挠度误差以 $O(h^{p+1})$ 下降。
也就是说，一次要素（$p=1$）网格尺寸减半时，挠度误差约降为1/4；二次要素（$p=2$）则约降为1/8。这个理论收敛率是否在实际网格收敛测试中复现，是代码验证（Code Verification）的核心。
实务上要用3种以上的网格密度求解，在对数-对数坐标上作图，看直线斜率是否接近理论值 $p+1$。如果偏离太大，说明有bug、边界条件错误或存在奇异点。
简支梁等分布荷载的实务应用
分析步骤
实际在求解器里跑这个基准问题，步骤是什么？
标准流程如下。
第1步：创建几何
$L \times b \times h = 1000 \times 100 \times 20$ mm 的直方体。可以利用对称性只建一半（$L/2$）。
第2步：生成网格
先从粗网格开始（跨径方向10分割，厚度方向2分割），逐步细化。收敛性测试建议跨径方向 10→20→40→80 分四个步骤。
第3步：定义材料
线性弹性。$E = 200{,}000$ MPa，$\nu = 0.3$。
第4步：设置边界条件
左端：$u_y = 0$（竖直方向约束）+ $u_z = 0$（面外位移约束）
右端：$u_y = 0$ + $u_z = 0$
其中一端加上 $u_x = 0$（约束轴向刚体位移）
第5步：施加荷载
上面施加面压 $p = w/b = 1.0/100 = 0.01$ MPa，或者施加线荷重 $w = 1.0$ N/mm。
第6步：求解与后处理
获取中央截面的最大挠度和最大应力，与理论解比对。
边界条件注意事项
边界条件设置有什么容易踩的坑？
常见的错误有3个。
1. 过度约束
「简支」本来就是两端各一个支点，但如果两端都设 $u_x = 0$，轴向伸长就被完全锁定了。可能出现热膨胀式的伪应力。应该只在一端设 $u_x = 0$。
2. 点约束引发应力集中
把支点约束在单个节点上，会在那点产生应力奇异。正确做法是沿支点所在的边或面约束一组节点，或用多点约束（MPC）定义支线。
3. 忘记面外方向约束
3D模型需要 $u_z = 0$ 来防止面外刚体转动。忘记的话刚性矩阵会奇异，求解器报错。
质量检查清单
分析完成后，怎样才能说「合格」？
基准测试合格的检查清单。
| 项目 | 合格条件 |
|---|---|
| 中央挠度 $\delta_{\max}$ | 与理论解误差 ±1% 以内 |
| 最大弯曲应力 $\sigma_{\max}$ | 与理论解误差 ±2% 以内 |
| 支点反力之和 | 等于 $wL$（力平衡验证） |
| 网格收敛性 | 两级以上细化显示单调收敛 |
| 收敛率斜率 | 对数-对数图的斜率在理论值 $p+1$ 的 ±0.2 内 |
| 沙漏能量比（低阶积分时） | 占总内部能的1%以下 |
| 对称性验证 | 关于中心左右对称的挠度分布 |
简支梁等分布荷载的软件比较
求解器各要素选择
各家求解器应该选什么要素？
主要求解器的推荐要素和设置整理如下。
| 求解器 | 推荐要素 | 自锁对策 | 输入卡/设置 |
|---|---|---|---|
| Abaqus | C3D20R（二次六面体低阶积分） | 内置 | *ELEMENT, TYPE=C3D20R |
| Abaqus | C3D8R（一次六面体低阶积分） | 沙漏控制 | *SECTION CONTROLS, HOURGLASS=ENHANCED |
| Ansys Mechanical | SOLID186（二次六面体） | 内置 | 默认设置即可 |
| Ansys Mechanical | SOLID185（一次六面体） | EAS法 | KEYOPT(2)=2 (Enhanced Strain) |
| MSC Nastran | CHEXA（20节点） | 内置 | PSOLID, FCTN=SMECH |
| CalculiX | C3D20R | 内置 | Abaqus兼容输入 |
| Code_Aster | HEXA20 | 内置 | MODELISATION='3D' |
交叉验证结果
求解器之间结果差异大吗？
采用二次六面体（HEX20同等品）的多求解器交叉验证结果。
| 求解器 | 要素 | $\delta_{\max}$ [mm] | 误差 [%] | $\sigma_{\max}$ [MPa] |
|---|---|---|---|---|
| Abaqus 2024 | C3D20R | 0.9762 | 0.04 | 18.73 |
| Ansys 2024R2 | SOLID186 | 0.9763 | 0.03 | 18.74 |
| MSC Nastran 2024 | CHEXA-20 | 0.9761 | 0.05 | 18.72 |
| CalculiX 2.22 | C3D20R | 0.9762 | 0.04 | 18.73 |
| 理论解 | — | 0.9766 | — | 18.75 |
全部求解器误差控制在0.05%以内。只要用二次要素，商用或开源都能得到相同精度。如果求解器间结果有较大偏差，就要检查要素定义或边界条件是否有问题。
简支梁等分布荷载的先端研究
非线性扩展
如果把这个问题拓展到非线性，会怎样？
两个扩展方向。
几何非线性（大变形）
当挠度接近或超过梁高 $h$ 时，曲率变化不能忽视，Euler-Bernoulli的线性假设破裂。此时需要von Karman型的大变形应变-位移关系。
当 $\delta/h > 0.5$ 时，非线性效应会超过5%，需要切换到这个定式。
材料非线性（弹塑性）
应力超过屈服应力时，需要von Mises屈服条件加上等向或运动硬化模型。在线性基准问题验证通过后，可进一步研究屈服载荷 $w_y$ 和极限载荷 $w_c$。
Richardson外推
没有解析解的问题怎样判断FEM结果的准确性？
用Richardson外推法。从多个网格密度的结果推算 $h \to 0$ 的极限值。
$f_1, f_2$ 是按细化顺序的网格解，$r$ 是网格细分比（通常为2），$p$ 是观测到的收敛次数。这个问题中有解析解，可以用Richardson外推值与解析解的一致性来验证该方法的有效性。掌握这个技巧后，实务中遇到没有解析解的复杂问题时就能用上了。
简支梁等分布荷载的故障排除
理论值不符情况
结果与理论值不符时，从何处检查？
排查的优先顺序。
1. 挠度过小（刚性过高）
- 剪切自锁 → 检查要素类型。完全积分的HEX8/TET4需要改成低阶积分，或升级到二次要素
- 过度约束 → 确认两端是否都被约束 $u_x = 0$
- 材料常数错误 → 确认 $E$ 的单位（GPa 和 MPa 的混淆是常见错误）
2. 挠度过大
- 荷载单位错误 → 确认 $w$ 是 N/mm 还是 N/m
- 尺寸输入错误 → $h$ 应该是 mm，如果误输为 m 会让 $I$ 偏差12个数量级
- 沙漏控制失效 → 低阶积分时，沙漏控制是否被正确设置
3. 只有应力值偏差
- 应力输出位置 → 积分点和节点外推值可能有差异。理论值比对时积分点最可靠
- 泊松比影响 → 3D固体中泊松比会产生微小的泊松效应（正常现象）
应力评估位置
应力应该在截面的哪里看？上面？下面？积分点还是节点？
理论上 $\sigma_{\max}$ 出现在中央截面（$x = L/2$）的最外缘（上面或下面）。FEM中有些细节要注意。
积分点（Gauss点）应力是最直接的计算结果，不含外推或平均化的误差。但积分点位于要素内部，不在表面。
节点外推值
实务上，网格足够细时，积分点值与节点外推值的差在1%以内。差异较大说明网格不够细。在基准问题中定量验证这个差异，就能培养在实际设计分析中判断应力值可信度的能力。 太有帮助了。简单问题仔细做，复杂问题就有底气了。 完全同意。简支梁基准就像「FEM的体温计」。掌握了要素特性、自锁现象、收敛规律，以后遇到任何问题都能问「这个结果靠不靠谱？」V&V的起点就是这个简单问题，强烈建议自己动手跑一遍，真实体验自锁和收敛。
有帮助
更详细
错误