3. 接触非线性赫兹接触分析
问题描述与建模思路
两个弹性球体在法向压力作用下形成局部接触区，接触半径远小于球体尺寸，属于典型的接触非线性问题。Abaqus/Standard 中采用接触对算法（主面-从面）处理该问题，接触压力分布与赫兹理论解对比是验证接触设置正确性的核心手段。
本例取两个相同钢球，半径 R=50 mm，弹性模量 E=210 GPa，泊松比 ν=0.3，法向载荷 P=1000 N。赫兹理论给出：接触半径 a=(3PR/4E*)^(1/3)，最大接触压力 p₀=3P/(2πa²)，其中 E*=E/(2(1-ν²))。
赫兹理论解计算
先算出理论值作为对比基准：
E* = E / (2*(1-ν²)) = 210e3 / (2*(1-0.09)) = 115384.6 MPa
a = (3*P*R / (4*E*))^(1/3)
= (3*1000*50 / (4*115384.6))^(1/3)
= (0.3250)^(1/3) ≈ 0.6875 mm
p₀ = 3*P / (2*π*a²)
= 3000 / (2*3.1416*0.4727) ≈ 1010.5 MPa理论接触半径约 0.6875 mm，最大接触压力约 1010.5 MPa，接触压力沿半径呈半椭球分布 p(r)=p₀√(1-(r/a)²)。
Abaqus 建模关键设置
- 几何：两个 1/4 球体（利用对称性），或完整球体取 1/8 模型；本例用 1/4 二维轴对称模型简化，半径 50 mm。
- 材料：线弹性，E=210000 MPa，ν=0.3。
- 分析步：Static, General，打开 NLGEOM=ON（接触本身即非线性）。
- 接触：Surface-to-surface contact，主面选下球上表面，从面选上球下表面；接触属性法向 Hard Contact，切向 Frictionless。
- 边界：下球底部固定，上球顶部施加集中力 P=1000 N（或位移控制）。
- 网格：接触区附近加密，单元尺寸约 0.05 mm，采用 CAX4I 或 C3D8I。
接触对算法中主从面的选择原则：刚度大、网格粗的一侧作主面。本例两球对称，可任选，但建议从面网格更细。
接触压力提取与对比
提交计算后，在 Job 结果中提取接触面法向压力 CPRESS，沿接触半径方向绘制曲线，与赫兹理论半椭球分布对比。判断标准：接触半径误差 <5%，最大压力误差 <5%，压力分布形状吻合。
结果判读与常见问题
| 检查项 | 理论值 | 可接受范围 | 常见偏差原因 |
|---|---|---|---|
| 接触半径 a | 0.6875 mm | 0.65~0.72 mm | 网格过粗、接触刚度不足 |
| 最大压力 p₀ | 1010.5 MPa | 960~1060 MPa | 接触面选择不当、单元类型 |
| 压力分布形状 | 半椭球 | 边缘平滑过渡 | 从面网格畸变、过约束 |
若接触压力出现振荡或尖峰，优先检查：①接触区网格是否足够细（建议单元尺寸≤a/10）；②是否使用了增强沙漏控制；③主从面是否发生穿透。可尝试调整接触刚度比例因子或改用罚函数法。
可复现操作清单
- 创建 2D 轴对称部件：上球为 1/4 圆，下球为 1/4 圆，半径 50 mm。
- 材料：E=210000 MPa，ν=0.3，截面赋给两个部件。
- 装配：两球顶点沿 Y 轴对齐，初始间隙 0（或微小过盈）。
- 分析步：Static, General，NLGEOM=ON，初始增量 0.01，最大 0.1。
- 接触：Surface-to-surface，主面下球上边，从面上球下边，法向 Hard，切向 Frictionless。
- 边界：下球底边 U2=0，上球顶边施加集中力 CF2=-1000 N。
- 网格：接触区种子 0.05 mm，其余 0.5 mm，单元 CAX4I。
- 提交 Job，提取 CPRESS 沿路径分布，与理论曲线对比。