布西内斯克问题（半无限弹性体的点荷载）

这是验证3D实体单元非常有用的经典问题。布辛涅斯克在1885年推导出了弹性半无限体表面施加集中荷载 $P$ 时的应力和位移精确解。因为它包含了3D场的$1/r$奇异性，所以非常适合评估单元精度和奇点附近的网格设计。

主要有两个。一个是用于FEA的3D实体单元的代码验证。另一个是赫兹接触分析的前期验证。因为赫兹接触理论是作为布辛涅斯克解的叠加（面压分布的卷积积分）推导出来的，所以用无法解决布辛涅斯克问题的求解器进行接触分析是危险的。

以荷载点为原点、荷载方向为$z$轴正方向，用圆柱坐标$(r, z)$表示。设$R = \sqrt{r^2 + z^2}$，位移场为

这里 $G = E/[2(1+\nu)]$ 是剪切弹性模量。应力场中最重要的$z$轴上的竖向应力为

这是个核心问题。FEA的离散化无法精确再现奇点。因此，需要在距离荷载点足够远的位置与理论值进行比较。评估点以 $z > 5h_{elem}$（$h_{elem}$是荷载点附近的单元尺寸）为基准。荷载点处的应力值因网格依赖而无意义，故不作为评估对象。

取足够大的模型尺寸即可。设关注区域的最大距离为 $L_{eval}$，则模型的外形应满足 $R_{model} \geq 10 L_{eval}$，深度 $D_{model} \geq 10 L_{eval}$。远场的边界条件标准做法是：底面固定$z$方向，侧面施加径向滚轮约束。

利用轴对称性，使用CAX8单元（Abaqus）或2D轴对称单元求解在计算上更高效。若需要3D验证，则使用1/4对称的3D模型。

在Abaqus中，在外缘布置CINAX4（轴对称无限单元）或CIN3D8（3D无限单元），可以将模型尺寸缩小到 $R_{model} \geq 3 L_{eval}$ 左右。计算成本会大幅下降，因此在进行参数化研究时推荐使用无限单元。

标准参数: $P = 1$ MN、$E = 200$ GPa、$\nu = 0.3$。评估点为荷载点正下方 $z = 0.1$ m。

在5mm网格尺寸下可获得1%以内的精度。用三个级别的网格计算GCI，结果约为0.5%，可以确认已充分收敛。

必须采用能应对$1/r^2$应力梯度的偏置网格。以荷载点为中心的放射状网格（蜘蛛网网格）是最优的，按几何级数增大单元尺寸。

• 荷载点附近: $h_{min} = L_{eval}/100$ 左右
• 外缘: $h_{max} = L_{eval}/2$ 左右
• 偏置比: 1.5〜2.0的几何级数

轴对称情况下，在荷载点（$r=0$轴上）布置一列退化的三角形单元，然后从该处向外放射状展开四边形单元。

精度可以达到，但达到同等精度所需的单元数是HEX20的3〜5倍。出于验证目的，推荐使用映射网格（六面体）。不过，以此问题进行面向复杂形状应用的自动四面体网格精度验证也很有意义。

以荷载点正下方$z = 0.1$ m处的$\sigma_{zz}$为指标，进行三个级别的计算。

观测收敛阶数: $p = \ln|(44.2-46.9)/(46.9-47.5)| / \ln 2 \approx 2.5$

GCI_fine: $1.25 \times |(-46.9-(-47.5))/(-47.5)| / (2^{2.5}-1) \approx 0.3\%$

在奇点附近，尚未完全进入渐近区域时可能会出现超收敛现象。应添加第四级网格（$h_{min} = 2.5$ mm）以确认$p$是否稳定。若能确认在渐近区域$p$稳定，则GCI的可靠性就得到了保证。

Abaqus: CAX8R是标准选择。用*CLOAD在轴上节点施加荷载。底面的无限单元用CINAX4。后处理可用Python脚本通过ODB API提取应力。

Nastran: 可用CTRIAX6进行轴对称分析但不好用。3D模型中使用CHEXA-20更实用。用SOL 101、GPSTRESS输出获取任意点的应力张量。

CalculiX: 使用与Abaqus兼容的C3D20。没有轴对称单元，所以用3D模型求解。用ccx执行，用ParaView可视化.frd文件。

Code_Aster: 用MODELISATION='AXIS'设置轴对称。使用相当于QUAD8的单元。用SALOME生成网格，在命令文件中设置边界条件。

在相同的网格密度和单元阶次下，位移可达到4〜5位有效数字一致，应力可达到3〜4位有效数字一致。差异主要源于节点外推算法的不同。若比较积分点处的值，则一致度会进一步提高。