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

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

主要有两个方面。一是FEA的3D固体单元的代码验证。另一个是Hertz接触分析的前期验证。Hertz接触理论是从Boussinesq解的叠加（面压分布的卷积积分）推导得出的，因此在不能通过Boussinesq问题的求解器上进行接触分析是危险的。

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

这是个核心问题。FEA的离散化无法精确再现奇异点。因此应在远离荷载点的位置进行与理论值的比较。评估点应以 $z > 5h_{elem}$（$h_{elem}$为荷载点附近的单元尺寸）为目标。荷载点处的应力值与网格有关，没有意义，因此不在评估范围内。

采用足够大的模型尺寸就可以了。如果关注领域的最大距离为 $L_{eval}$，模型外形应为 $R_{model} \geq 10 L_{eval}$，深度 $D_{model} \geq 10 L_{eval}$。远方边界条件的标准做法是底面$z$方向固定，侧面为径向滚轮。

可以的。使用Abaqus的CINAX4（轴对称无限单元）或CIN3D8（3D无限单元）放在外边缘，可以将模型尺寸缩减到$R_{model} \geq 3 L_{eval}$左右。计算成本大幅降低，所以在参数化研究中建议采用无限单元。

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

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

精度可以，但获得相同精度所需的单元数是HEX20的3～5倍。对于验证目的，建议使用映射网格（六面体）。不过，从为复杂形状应用考虑出发，对自动四面体网格的精度进行此问题验证也很有意义。

以荷载点正下方$z = 0.1$ m的$\sigma_{zz}$为指标进行3个水平的计算。

在奇异点附近，在进入完全渐近区域之前可能出现超收敛。应添加第4个网格水平（$h_{min} = 2.5$ mm）来验证$p$是否稳定。一旦在渐近区域内$p$稳定，就可保证GCI的可靠性。

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

在相同网格密度和单元次数下，位移的4～5位有效数字相同，应力的3～4位相同。差异主要来自节点外推算法的差异。如果在积分点值处进行比较，一致度会进一步提高。

按以下步骤进行。

评估点($z = 0.1$ m)处$\sigma_{zz}$和$u_z$与理论值的偏差在1%以内。以"细"网格的值为基准，GCI < 5%。此基准与悬臂梁基准处于同一精度水平，具有一致性。

Hertz接触理论是Boussinesq解对面压分布的卷积积分推导得出的。球与平面接触时，接触面处的半椭圆形面压分布

完全正确。V&V的基本原则是"从简单到复杂"。Boussinesq（线性弹性、轴对称）→Hertz（接触但有解析解）→一般接触（非线性、有摩擦）的推进。在各阶段确认与参考解的一致，再以此作为下一阶段的信任基础。

在地基设计初期作为第一近似被使用。计算地中增加应力，推估压密沉降量。具体采用Newmark图表或Fadum积分表的方法来求取矩形分布荷载下的应力增量，这种方法已广泛应用。

当作用应力远小于地基屈服应力的地域（常时荷载下的地基）时，弹性近似就能获得实用精度。如需处理大变形或液化，则需进入弹塑性模型或有效应力分析。但即使这样，Boussinesq解在初期应力状态推估或数值级估计上仍可使用。

轴对称CAX8R模型的要点如下。

*OUTPUT, HISTORY, VARIABLE=PRESELECT 输出特定节点的时刻历位移。应力用*EL PRINT 获取单元积分点值，或用Python ODB API提取坐标指定的场输出。probeValues 方法可获取任意点的内插值，非常便于理论值自动对比。

CTRIAX6是6节点三角形，四边形的情况下用TRAX3或TRAX6，但不如Abaqus方便。实务中常用3D 1/4对称模型配合CHEXA-20求解。

Nastran的应力张量输出顺序为 $\sigma_{xx}, \sigma_{yy}, \sigma_{zz}, \sigma_{xy}, \sigma_{yz}, \sigma_{xz}$，而Abaqus为 $\sigma_{11}, \sigma_{22}, \sigma_{33}, \sigma_{12}, \sigma_{13}, \sigma_{23}$。特别是剪切分量的顺序不同（Nastran为 $xy, yz, xz$，Abaqus为 $12, 13, 23$），对比时一定不要弄混。

CalculiX与Abaqus的.inp文件格式基本兼容，几乎可直接使用。用C3D20单元进行3D 1/4模型求解。需要注意的是，CalculiX的C3D20积分点配置与Abaqus略有不同。结果精确到3～4位，但要求更高精度时应在源代码（src/e_c3d.f）中确认积分方案细节。

源代码公开，第三方能直接验证单元定式推导算法。实际上CalculiX的C3D20刚度矩阵计算在e_c3d.f内的Gauss积分中明确呈现。商业求解器则没有这种透明性。在学术论文中基于开源软件进行验证，其结果可作为商业求解器的独立验证引用，这是既定方法论。

Burmister（1945）推导了双层弹性体解。其中包含Hankel变换的积分形式，需数值积分但作为参考解存在。铺装设计（上层沥青、下层路基）的FEA验证中使用。

完全粘结条件下共享节点足够，但若考虑滑动界面则需接触条件。完全粘结的验证以Burmister解为参考；滑动界面的验证需用层间无摩擦的Burmister解。分阶段验证尤为重要。

Lamb问题（1904）正是这种情况。半无限体表面施加阶跃荷载时的过渡性弹性波传播，有P波、S波、Rayleigh波到达时刻与振幅的理论解。

每波长的节点数是关键。二次单元1波长5～6个单元，一次单元10～12个单元为目标。S波波长 $\lambda_s = V_s / f$ 决定所需单元尺寸。需对不同网格密度的波形进行比较，确认数值色散（周频依赖的相速度偏差）在允许范围内。

弹塑性无精确解，以Boussinesq弹性解作为杨氏模量参数化研究的基准。低荷载阶段弹塑性解与弹性解一致，随荷载增加偏离扩大（压缩屈服→塑性域扩展→承载力极限）。确认此偏离规律物理合理即可。

正是如此。以简单问题确立基础信任，逐步复杂化。Boussinesq → Hertz接触 → 多层地基 → 弹塑性 → 动力问题的系统化V&V金字塔中，Boussinesq居其底层。这种系统方法正是ASME V&V 10和NAFEMS QSS推荐的方法论。

固定边界的反射导致应力被高估。底面固定过浅时，荷载点正下方的$\sigma_{zz}$会比理论值大20%以上。现象为改变模型尺寸时结果波动，则应怀疑此影响。

无限单元的衰减函数可能与问题特性不匹配。Abaqus无限单元假设指数衰减，但Boussinesq解是$1/r$衰减，严格上不完全吻合。实用上精度足够，但若对残差敏感，加大有限单元区域（向无限单元边界远离）较为安全。

这是正常现象。荷载点处的应力随网格细化无限增长。这反映了物理的奇异性，不是FEA缺陷。

用低阶单元（TET4、HEX8）的应力等值线图会产生锯齿形不连续。这是因为C0连续单元在元素边界处应力（歪变导数）不连续。

Zienkiewicz和Zhu（1992）提出的事后误差估计法。各节点的斑块内积分点应力用最小二乘法高精度恢复为节点应力。恢复前后的差是局部离散化误差估计。商用求解器中多作为选项实现。Abaqus中可用*ERROR INDICATOR。