平行圆柱体的赫兹接触计算与ANSYS实现

CAE_LJX 2021年3月30日浏览：2442 评论：28 收藏：15

关注公众号：“CAE之道”，享受专属答疑服务，精彩文章不错过。

1882年，年仅25岁的德国天才物理学家赫兹发表了关于接触力学的著名文章《关于弹性固体的接触（On the contact of elastic solids）》，系统地阐述了两物体之间接触面上所传递的压力分布，以及它所引起的垂直于接触面的弹性位移在接触区内、外的关系。另外，赫兹在这篇论文中提出了有关弹性体接触的理论公式——赫兹公式。

赫兹

海因里希·鲁道夫·赫兹是德国天才物理学家，毕业于柏林大学，师承古斯塔夫·基尔霍夫（德国物理学家，在电路、光谱学的基本原理有重要贡献）和赫尔曼·范·亥姆霍兹（能量守恒定律的创立者）。

赫兹在波动方程、光电效应和接触力学等领域都做出了杰出的贡献，尤其是他于1888年首次通过实验证明了电磁波的存在，成为了近代科学史上的一座里程碑，开创了无线电电子科学技术的新纪元。

正当人们对他寄以更大期望时，他却于1894年元旦因血中毒逝世，年仅36岁。为了纪念他的功绩，人们用他的名字来命名各种波动频率的单位—赫兹（Hz）。频率也是国际单位制中频率的单位，它是每秒中周期性变动重复次数的计量。

图为赫兹及夫人伊丽莎白

赫兹公式是研究疲劳、摩擦以及任何有接触体之间相互作用的基本公式。接触理论指出：接触表面上所承受的压应力是处处不同的，其分部呈半椭圆柱形。初始接触线处压应力最大，以此最大压应力代表两零件间接触受力后的应力。

平行圆柱体的赫兹接触计算与ANSYS实现的图2

平行圆柱体的赫兹接触计算与ANSYS实现的图3

赫兹公式也是基于一定的假设，其作出的假设如下：

用a表示接触区的有效尺寸，用ρ表示曲率半径，用R表示每个物体的有效半径，用l表示物体横向和深度两方面的有效尺寸，则赫兹理论中做出的假设可以简单表述成：

1. 表面都是连续的，并且是非协调的：a〈〈 ρ；

2. 接触尺寸远小于接触物体尺寸；

3. 小应变；

4. 每个接触物体都是线弹性的，服从胡克定律；

5. 接触物体间摩擦力为0。

为了对赫兹公式的计算结果和ANSYS的计算结果进行对比，我们选择以两横截面直径为100mm、b为100mm，泊松比为0.3、弹性模量为200Gpa的 长圆柱体为例，假设外载F=20kN，分别基于 赫兹公式和 ANSYS软件计算一下接触面面半宽和最大接触应力：

一、基于赫兹公式的计算：

为了计算方便，此处笔者将赫兹公式编制成了一个简单的Python小程序，代码及计算结果如下：

平行圆柱体的赫兹接触计算与ANSYS实现的图4

根据计算结果我们发现，该问题中两物体的接触面半宽为0.2407mm，远小于接触物体的结构尺寸，因此 符合赫兹公式的假设。

二、基于ANSYS软件的计算：

使用ANSYS求解该问题时，我们从以下几个方面入手：

1. 确定分析类型：根据例题所示结构，确定分析类型为 静力学分析；

2. 确定单元类型： 两长圆柱体的分析计算，为了降低计算量，可使用1/4的平面应变模型计算（具体选用规则请看本公众号《ANSYS与材料力学之轴向拉伸和压缩（二）》）。

Step1

平面模型建模

在保证精度的前提下，为了降低计算量，本次分析使用1/4的平面模型，建模方式如下：

平行圆柱体的赫兹接触计算与ANSYS实现的图6

Step2

建立分析模块

打开Workbench，选择Static Structural模块，并传入上一步建立的几何模型。

平行圆柱体的赫兹接触计算与ANSYS实现的图7

Step3

平面几何设置

在Project Schematic中的空白处点击右键，选择Properties，打开Properties of Project Schematic。单击项目中的A2（Geometry）栏，在Propertiesof Project Schematic A2: Geometry中将AnalysisType切换为2D。

平行圆柱体的赫兹接触计算与ANSYS实现的图8

Step4

材料设置

为了保证与理论计算输入参数的一致性，本次分析使用默认的 结构钢材料（弹性模量为200000MPa，泊松比为0.3）。

平行圆柱体的赫兹接触计算与ANSYS实现的图9

Step5

平面分析设置

在结构树中点击Geometry，将Details of Geometry中的2D Behavior切换成Plane Strain（平面应变）。

Step6

轴对称设置

由于我们建模的时候使用了1/4的模型，所以我们需要在对称轴处设置对称关系。

建立一个轴对称区域，选择上下两个平面模型的对称轴，并将对称法向设置为x轴。

平行圆柱体的赫兹接触计算与ANSYS实现的图10

Step7

接触设置

1. 接触边和目标边分别选择两段圆弧边；

2. 接触类型选择Frictionless（无摩擦接触）；

3. 接触算法选择Aaugmented Lagrange（增广拉格朗日算法）；

4. 法向刚度因子设置为20。

Step7

网格划分

1. 全局网格尺寸设置为2mm，一阶单元；

2. 使用影响球对接触位置进行局部加密。此时我们需要在要局部加密的位置建立一个局部坐标系。

3. 由于我们计算出的接触面半宽为0.2407mm，因此我们把影响球半径设置为0.5mm，网格尺寸设置为0.01mm，其余保持默认。

平行圆柱体的赫兹接触计算与ANSYS实现的图14

Step9

求解设置

1. 打开自动时间步；

2. 时间步类型选择 Substeps（子步）；

3. 初始子步设置为 50；

最小子步设置为 50 ；

最大子步设置为100。

平行圆柱体的赫兹接触计算与ANSYS实现的图15

Step10

约束设置

设置底边为 固定约束。

平行圆柱体的赫兹接触计算与ANSYS实现的图16

Step11

载荷设置

理论计算时载荷为20kN，我们现在使用的是平面应变模型，因此需要对载荷进行一定的换算。

此处笔者使用的是 Pressure加载，具体的换算方法是：

Pressure = F /（圆柱直径*圆柱长度）

= 20000 /（100*100）

= 2MPa

Step12

求解及后处理

在结果中插入 接触工具Contact Tool，提取 接触压力Pressure（赫兹公式中的接触应力），可以看出，计算结果为 528.56MPa，与赫兹公式解出的 528.9923MPa 几乎完全一致。

平行圆柱体的赫兹接触计算与ANSYS实现的图18

对于接触面半宽的提取，我们可以在ANSYS APDL中实现，具体步骤如下：

Step13

建立Workbench与APDL的连接

平行圆柱体的赫兹接触计算与ANSYS实现的图19

注意：如果在Workbench的分析中使用了上述 Step6中的轴对称条件，在建立连接时，软件会报错，错误类型如下图所示。解决办法是去掉 Step6， 设置约束时在对称边上建立无摩擦约束。

平行圆柱体的赫兹接触计算与ANSYS实现的图20

Step14

查看接触压力分布

平行圆柱体的赫兹接触计算与ANSYS实现的图21

在通用后处理的Results Viewer中，显示接触压力结果。我们发现接触部分的接触压力的确呈半椭圆分布。

Step15

查看接触面半宽

笔者查看接触面半宽的方式是：使用 DISTNP（N1,N2）函数 测量最大接触压力节点和最小接触压力节点的距离。通过接触压力的分布图，我们找到最大接触压力节点为3141节点，最小接触压力节点为3227节点。我们在命令行中输入b = DISTND（3141,3227），命令输出窗口即显示b的值为 0.2488。

对比使用赫兹公式计算出的接触面半宽0.2407mm，ANSYS计算的接触面半宽0.2488mm， 误差为3.4%，也是可以接受的。

平行圆柱体的赫兹接触计算与ANSYS实现的图22