Abaqus 带孔薄板拉伸：应力集中与网格精度

1 横向均匀拉伸矩形带孔薄板问题

  如下图所示，一个厚度均匀的矩形薄板，左右两边受均匀拉伸载荷作用。已知薄板尺寸为：长 $a=1m$，宽 $b=0.6m$，厚 $t=0.01m$，在薄板中心位置挖掉一个半径为 $r=0.1m$ 的圆孔。薄板左右两边均受 $q=100MPa$ 的均匀载荷作用。薄板的弹性模量 $E=200GPa$，泊松比 $\mu=0.3$，忽略自重。试使用不同的单元类型及不同的单元数量（至少选择 $3$ 组不同的情况）对薄板进行有限元分析（应力、应变及位移等），并与解析解对比计算精度。

  本文中使用的Abaqus版本是2025。

2 Abaqus模型构建

  Abaqus源文件点此处下载。

3 网格划分与数值模拟云图

  共计算三组：六面体网格（全局尺寸 0.05）、四面体网格（全局尺寸 0.05）、四面体网格（全局尺寸 0.005）。每组给出网格划分与应力、应变、位移云图。

3.1 六面体网格，全局尺寸0.05

3.2 四面体网格，全局尺寸0.05

3.3 四面体网格，全局尺寸0.005

4 数值模拟结果与解析结果比较

  先进行理论计算。由基尔斯（Kirsch）解答可知 $\sigma_r$、$\sigma_\theta$、$\tau_{r\theta}$ 的表达式，其中 $a$ 为圆孔半径：

$$\begin{aligned} \begin{cases} \sigma_r=\dfrac{q}{2}\left(1-\dfrac{a^2}{r^2}\right)+\dfrac{q}{2}\left(1-\dfrac{4a^2}{r^2}+\dfrac{3a^4}{r^4}\right)\cos2\theta\\ \sigma_\theta=\dfrac{q}{2}\left(1+\dfrac{a^2}{r^2}\right)-\dfrac{q}{2}\left(1+\dfrac{3a^4}{r^4}\right)\cos2\theta\\ \tau_{r\theta}=-\dfrac{q}{2}\left(1+\dfrac{2a^2}{r^2}-\dfrac{3a^4}{r^4}\right)\sin2\theta \end{cases} \end{aligned}$$

  孔边应力集中点取 $r=a$、$\theta=\dfrac{\pi}{2}$、$q=100MPa$，代入得：

$$\begin{aligned} \begin{cases} \sigma_r=0\\ \sigma_\theta=3q=300MPa\\ \tau_{r\theta}=0 \end{cases} \end{aligned}$$

  该点处于与拉伸方向垂直的位置，化到直角坐标系即：

$$\begin{aligned} \begin{cases} \sigma_x=300MPa\\ \sigma_y=0\\ \tau_{xy}=0 \end{cases} \end{aligned}$$

  与第3章结果进行比较，见下表。

5 总结与讨论

  本次以带圆孔矩形薄板为对象，用 Abaqus 计算三组不同网格，并与基尔斯解比较。结果汇总如下表。

  网格类型的影响：在相同网格尺寸（$0.05m$）下，六面体网格的 $\sigma_x$ 相对误差为 $7.34\%$，低于四面体网格的 $8.60\%$；$\tau_{xy}$ 更是从 $52.20MPa$ 降到 $12.05MPa$（理论值为 0）。这与两类单元的固有性质一致——六面体单元形函数插值精度更高，且在规则几何上更容易生成形态良好的单元。

  网格尺寸的影响：把四面体网格从 $0.05m$ 细化到 $0.005m$（单元尺寸缩小到 1/10）后，三个分量的表现并不一致。$\sigma_y$ 的误差由 $24.06MPa$ 降到 $0.34MPa$，$\tau_{xy}$ 由 $52.20MPa$ 降到 $4.64MPa$——两者都大幅改善。$\sigma_x$ 则从 $325.81MPa$ 升到 $355.212MPa$，相对误差由 $8.60\%$ 变为 $18.40\%$。

  这三个分量的差异，恰好可以用同一个机制解释清楚。$\sigma_y$ 与 $\tau_{xy}$ 的理论值都是 0，它们属于“孔边的非峰值分量”，本来就应该趋于 0；网格细化使单元内的应力变化减小，读数自然更贴近 0，所以误差大幅下降。$\sigma_x$ 的理论值则是孔边峰值 $300MPa$——孔边是应力集中点，$\sigma_x$ 从孔边的峰值向外迅速衰减，而有限元给出的是单元内的近似（本例还采用了 75% 的节点应力平均，进一步平滑了峰值）。单元越粗，单元覆盖的范围内应力落差越大，平均值就越低于真实峰值，因此粗网格低估 $\sigma_x$；网格细化后单元覆盖范围缩小，平均值趋近峰值，所以 $\sigma_x$ 随细化升高。

  换句话说：$\sigma_y$、$\tau_{xy}$ 的读数向 0 靠拢，$\sigma_x$ 的读数向 $300MPa$ 的峰值方向抬升，两者的变化方向不同，但都是各自向真值收敛的正确表现。$\sigma_x$ 升高不是错误，而正是细化在起作用。

  至于三组 $\sigma_x$ 与理论值仍有差距（$277.98$、$325.81$、$355.212MPa$ 对应 $300MPa$），原因是网格相对孔径仍不够细。本例单元尺寸 $0.005m$ 与孔径 $0.1m$ 之比为 $1/20$，而应力集中点的峰值需要远细于此的局部网格（通常要求孔边单元尺寸达到孔径的 $1/50$ 甚至更小）才能准确捕捉；再叠加 75% 节点平均带来的平滑，读数与理论峰值存在百分之十几的差距属正常范围。若继续在孔边做局部加密，$\sigma_x$ 会进一步趋近 $300MPa$ 并趋于稳定，这是应力集中问题的标准收敛过程。因此本文结论是：对应力集中问题，网格细化会提升峰值应力的求解精度，但要求孔边局部单元尺寸远小于孔径，全局网格加密的代价高而收益有限。

  位移分布规律：由位移云图与数值输出可见，$U1$（水平）位移随网格细化轻微增大（六面体 $0.05m$ 下约 $2.1\times10^{-3}m$，四面体 $0.05m$ 下 $2.3\times10^{-3}m$，四面体 $0.005m$ 下 $2.6\times10^{-3}m$），反映细化后对局部变形的捕捉更充分；$U2$（竖直）位移在三种方案下均在 $1.5\sim1.7\times10^{-3}m$ 之间，波动很小，说明它主要由整体受力特征决定；$U3$（厚度方向）位移在 $10^{-6}m$ 量级，可忽略——这从数值上印证了把该问题按平面应力简化的合理性。

  关于左右不对称现象：计算中发现位移云图的形状左右对称，但数值并不左右对称，多组重复计算均如此。一种可能是两侧均布载荷的施加顺序（先右后左）引入了整体的刚体位移分量；若确实如此，可用“变形前孔心与变形后孔心的水平距离”作为修正量，从各点水平位移中扣除。这一点尚未做独立的验证实验（例如把载荷改为一次施加、或改用位移边界），因此这里只作为现象记录，不作为结论。