简支梁均布荷载 — FEM验证基准

这是FEM最基础的验证问题。承受均布荷载 $w$ 的简支梁，根据欧拉-伯努利梁理论可以得到精确的挠度解。中央的最大挠度为

将理论解与数值解进行比较，以确认单元类型的精度，是这个基准测试的目的。特别重要的是，使用一阶单元（线性单元）会发生剪切自锁，导致梁表现得比实际更硬，挠度被低估。定量地把握这一现象，并验证减缩积分或B-bar法的改善效果，是实践中的关键点。

是的。正因为简单才好用。它被收录在NAFEMS的基准测试集中，也作为ASME V&V 10-2006中代码验证（确认代码数学正确性）的第一步被推荐。引入新求解器时，首先用这个问题确认与理论解在0.1%以内的一致，是业界的标准流程。例如，汽车制造商的构造解析部门，在每次求解器版本升级时，都会将这个基准测试作为回归测试来运行。

从欧拉-伯努利梁的控制方程开始吧。承受均布荷载 $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$，在中央为零。欧拉-伯努利理论假设截面保持平面且法线方向不变（忽略剪切变形），因此对于跨度/高度比 $L/h$ 较大的薄梁精度较高。

铁木辛柯梁理论额外考虑了剪切变形。中央挠度为

第二项是剪切变形引起的挠度增量。$\kappa$ 是剪切修正系数（矩形截面为 $5/6$），$G$ 是剪切弹性模量 $G = E/\{2(1+\nu)\}$。

若 $L/h \geq 20$，则欧拉-伯努利理论已足够。对于 $L/h < 10$ 的深梁，则需要铁木辛柯理论或3D实体单元。基准测试通常使用 $L/h = 50$ 左右的条件，在欧拉-伯努利理论严格成立的条件下进行验证。

标准的基准测试设置如下。

$L/h = 50$，因此欧拉-伯努利理论充分成立。截面惯性矩 $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$。

使用上述参数的比较结果如下。

HEX8完全积分出现了约17%的位移误差，TET4甚至高达30%。这就是剪切自锁的威力。另一方面，减缩积分的HEX8或HEX20误差在0.3%以内。清晰地展示这种差异，正是简支梁基准测试的意义所在。

要理解本质，可以考虑一阶单元形函数的局限性。

一阶六面体单元（HEX8）的位移场在各个方向都是线性的。纯弯曲时截面会旋转，因此上表面和下表面会产生相反方向的水平位移。理论上此时剪切应变 $\gamma_{xy} = 0$，但线性的位移场无法满足这个条件。

也就是说，产生了本不该存在的寄生剪切应变（parasitic shear）。对应于该应变的剪切应力会作为额外的刚度起作用，使得梁表现得比实际更硬，即挠度被低估。这就是剪切自锁。

没错。即使细化网格，单个单元的纵横比得到改善，情况会缓解，但收敛速度非常慢。厚度方向只有2个单元左右是远远不够的，需要大约8~16个单元。计算成本会变得巨大，因此在实践中使用自锁对策技术要高效得多。

主要有三种方法。

1. 减缩积分（Reduced Integration）
使用减缩积分（1×1×1=1点）代替完全积分（HEX8为2×2×2=8点）。通过减少积分点，可以排除寄生剪切应变的能量贡献。但代价是会产生沙漏模式（hourglass mode）这种零能量变形模式。即使像沙漏那样扭曲变形，也会被评估为零能量。

2. B-bar法（选择性减缩积分）
仅对体积变化（膨胀）分量使用减缩积分计算，而偏差分量则保持完全积分。对体积自锁和剪切自锁都有效。Abaqus的 C3D8R + hourglass control 或 Nastran的 CHEXA reduced 就是基于这个思路。

3. EAS法（增强假设应变）
在通常的位移场之外，假设单元内部存在额外的应变模式。引入内部自由度作为额外的应变参数，并通过静态凝聚使其对外部不可见。Ansys的SOLID185（keyopt(2)=2）就属于此类。

通过添加人工刚度（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$），将单元尺寸 $h$ 减半，位移误差大约变为 $1/4$；对于二阶单元（$p=2$），则大约变为 $1/8$。在实际的网格收敛测试中确认是否能重现这个理论收敛率，是代码验证的核心。

具体来说，使用三种以上不同的网格密度进行分析，在双对数图上确认斜率是否接近理论值 $p+1$。如果斜率与理论值偏差较大，则需要怀疑是否存在错误、边界条件设置不当或奇异点的影响。

基本流程如下。

步骤 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: 求解与后处理
比较中央截面的最大挠度和最大应力与理论解。

常见的陷阱有三个。

1. 过约束
明明是"简支"，却在两端都设置了 $u_x = 0$ 的情况。轴向的热膨胀之类的伸长被约束，会产生寄生的轴向力。应该只在一端设置 $u_x = 0$。

2. 点约束引起的折角
在实体单元中，如果仅用一个点（单一节点）约束支点，那里会产生应力集中。正确的做法是沿着支点侧的边或面约束节点群，或者用MPC（多点约束）定义支撑线。

3. 忘记约束面外方向
在3D模型中，为了防止 $z$ 方向的刚体旋转，需要在面外方向设置 $u_z = 0$。如果忘记设置，会导致刚度矩阵奇异，求解器报错。

这是基准测试合格的检查清单。