赫兹接触理论（球-平面）

因为对于两个弹性体的接触，接触半径、最大接触压力、接触位移等所有量都存在闭合形式的解析解。这是赫兹在1882年推导出来的。球与平面的接触是最基本的，是非线性接触分析代码验证的入门。

赫兹理论的巧妙之处就在于此。它包含了接触面积随载荷变化的几何非线性，但通过联立弹性半空间的布辛涅斯克解叠加和几何相容条件，可以得到闭合形式解。不过有前提条件：接触面相对于平面要足够小、在弹性限度内、无摩擦。

半径为 $R$ 的弹性球以载荷 $P$ 压向弹性平面的情况。定义等效弹性模量 $E^* = [(1-\nu_1^2)/E_1 + (1-\nu_2^2)/E_2]^{-1}$，等效半径 $R^* = R$（平面 $R_2 = \infty$），则有

接触半径: $a = \left(\frac{3PR^*}{4E^*}\right)^{1/3}$

最大接触压力: $p_0 = \frac{2aE^*}{\pi R^*} = \frac{1}{\pi}\left(\frac{6PE^{*2}}{R^{*2}}\right)^{1/3}$

趋近量: $\delta = \frac{a^2}{R^*} = \left(\frac{9P^2}{16R^*E^{*2}}\right)^{1/3}$

接触面内的压力呈半椭圆分布。

中心处最大为 $p_0$，接触边缘 $r = a$ 处为零。这个分布可以通过布辛涅斯克解的叠加来再现，接触面下的应力场也能精确计算。最大 von Mises 应力出现在表面下方约 $z \approx 0.48a$ 处，是滚动疲劳的起点。

钢球（$R = 10$ mm、$E = 200$ GPa、$\nu = 0.3$）对钢平面（相同材料）。$P = 100$ N。

$E^* = 200/(2 \times (1-0.09)) = 109.9$ GPa

$a = (3 \times 100 \times 0.01 / (4 \times 109.9 \times 10^9))^{1/3} = 0.0727$ mm

$p_0 = 3P/(2\pi a^2) = 9013$ MPa

$\delta = a^2/R = 0.000529$ mm

这是在接触面内布置了20个以上单元的结果。因为接触半径 $a$ 非常小，接触区域的单元尺寸需要设为 $a/10 \approx 7$ μm。远离区域可以粗化，所以必须使用偏置网格。

主要有3种。

对于赫兹验证，推荐使用拉格朗日乘子法或增广拉格朗日法。罚函数法会因为罚函数系数的取值而改变接触压力，不适合用于验证。

Abaqus 默认使用增广拉格朗日法（SURFACE INTERACTION + CONTACT PAIR）。Nastran 的 SOL 101 中 GLUE/SLIDE 接触基于罚函数。Ansys 默认使用增广拉格朗日法。在接触问题的验证与确认中，必须注明所使用的算法。

接触面的单元尺寸必须是接触半径 $a$ 的1/10以下。网格太粗会使接触面呈阶梯状，导致接触压力分布不准确。

推荐配置:

• 接触面附近: $h_{elem} = a/10$。使用 HEX20 结构网格
• 远离接触面的区域: 使用偏置网格粗化。比率约 2.0
• 球的1/4对称模型（两个对称面）
• 平面侧模型尺寸取半径 $10a$ 以上

按照 Abaqus 的惯例，凸面（球）设为从面，平面设为主面。从面的网格应该更细。Nastran 和 Ansys 也采用类似思路，将刚度更高的面设为主面。如果两侧网格密度差异很大，会影响接触压力的精度，因此最好在接触面处使两侧网格密度保持一致。

像赫兹问题这样光滑的接触通常能顺利收敛。如果不收敛：

1. 检查初始接触间隙。间隙过大会导致初始步无法检测到接触

2. 将载荷分成多个载荷步。启用 NLGEOM=ON 以考虑几何非线性

3. 临时启用接触稳定化（Abaqus 的 *CONTACT CONTROLS, STABILIZE）

4. 使用平滑选项（Abaqus 的 *SURFACE SMOOTHING）减轻接触面的不连续性

在赫兹接触中，球与平面的趋近伴随着接触面积的变化。这是几何非线性，因此必须设置 NLGEOM=ON。如果不开启，求解器会基于初始间隙状态进行线性近似，无法得到正确的接触面积。