平行圆柱体的赫兹接触计算与ANSYS实现
1882年,年仅25岁的德国天才物理学家赫兹发表了关于接触力学的著名文章《关于弹性固体的接触(On the contact of elastic solids)》,系统地阐述了两物体之间接触面上所传递的压力分布,以及它所引起的垂直于接触面的弹性位移在接触区内、外的关系。另外,赫兹在这篇论文中提出了有关弹性体接触的理论公式——赫兹公式。
海因里希·鲁道夫·赫兹是德国天才物理学家,毕业于柏林大学,师承古斯塔夫·基尔霍夫(德国物理学家,在电路、光谱学的基本原理有重要贡献)和赫尔曼·范·亥姆霍兹(能量守恒定律的创立者)。
赫兹在波动方程、光电效应和接触力学等领域都做出了杰出的贡献,尤其是他于1888年首次通过实验证明了电磁波的存在,成为了近代科学史上的一座里程碑,开创了无线电电子科学技术的新纪元。
正当人们对他寄以更大期望时,他却于1894年元旦因血中毒逝世,年仅36岁。为了纪念他的功绩,人们用他的名字来命名各种波动频率的单位—赫兹(Hz)。频率也是国际单位制中频率的单位,它是每秒中周期性变动重复次数的计量。图为赫兹及夫人伊丽莎白
赫兹公式是研究疲劳、摩擦以及任何有接触体之间相互作用的基本公式。接触理论指出:接触表面上所承受的压应力是处处不同的,其分部呈半椭圆柱形。初始接触线处压应力最大,以此最大压应力代表两零件间接触受力后的应力。
赫兹公式也是基于一定的假设,其作出的假设如下:
用a表示接触区的有效尺寸,用ρ表示曲率半径,用R表示每个物体的有效半径,用l表示物体横向和深度两方面的有效尺寸,则赫兹理论中做出的假设可以简单表述成:
1. 表面都是连续的,并且是非协调的:a〈〈 ρ;为了对赫兹公式的计算结果和ANSYS的计算结果进行对比,我们选择以两横截面直径为100mm、b为100mm,泊松比为0.3、弹性模量为200Gpa的长圆柱体为例,假设外载F=20kN,分别基于赫兹公式和ANSYS软件计算一下接触面面半宽和最大接触应力:为了计算方便,此处笔者将赫兹公式编制成了一个简单的Python小程序,代码及计算结果如下:
根据计算结果我们发现,该问题中两物体的接触面半宽为0.2407mm,远小于接触物体的结构尺寸,因此符合赫兹公式的假设。二、基于ANSYS软件的计算:
使用ANSYS求解该问题时,我们从以下几个方面入手:
1. 确定分析类型:根据例题所示结构,确定分析类型为静力学分析;在保证精度的前提下,为了降低计算量,本次分析使用1/4的平面模型,建模方式如下:
打开Workbench,选择Static Structural模块,并传入上一步建立的几何模型。
在Project Schematic中的空白处点击右键,选择Properties,打开Properties of Project Schematic。单击项目中的A2(Geometry)栏,在Propertiesof Project Schematic A2: Geometry中将AnalysisType切换为2D。
为了保证与理论计算输入参数的一致性,本次分析使用默认的结构钢材料(弹性模量为200000MPa,泊松比为0.3)。
在结构树中点击Geometry,将Details of Geometry中的2D Behavior切换成Plane Strain(平面应变)。
由于我们建模的时候使用了1/4的模型,所以我们需要在对称轴处设置对称关系。
建立一个轴对称区域,选择上下两个平面模型的对称轴,并将对称法向设置为x轴。
2. 接触类型选择Frictionless(无摩擦接触);
3. 接触算法选择Aaugmented Lagrange(增广拉格朗日算法);
4. 法向刚度因子设置为20。
2. 使用影响球对接触位置进行局部加密。此时我们需要在要局部加密的位置建立一个局部坐标系。3. 由于我们计算出的接触面半宽为0.2407mm,因此我们把影响球半径设置为0.5mm,网格尺寸设置为0.01mm,其余保持默认。
最大子步设置为100。
理论计算时载荷为20kN,我们现在使用的是平面应变模型,因此需要对载荷进行一定的换算。
此处笔者使用的是Pressure加载,具体的换算方法是:Pressure = F /(圆柱直径*圆柱长度) = 20000 /(100*100)
在结果中插入接触工具Contact Tool,提取接触压力Pressure(赫兹公式中的接触应力),可以看出,计算结果为528.56MPa,与赫兹公式解出的528.9923MPa几乎完全一致。
对于接触面半宽的提取,我们可以在ANSYS APDL中实现,具体步骤如下:
注意:如果在Workbench的分析中使用了上述Step6中的轴对称条件,在建立连接时,软件会报错,错误类型如下图所示。解决办法是去掉Step6,设置约束时在对称边上建立无摩擦约束。
在通用后处理的Results Viewer中,显示接触压力结果。我们发现接触部分的接触压力的确呈半椭圆分布。笔者查看接触面半宽的方式是:使用DISTNP(N1,N2)函数测量最大接触压力节点和最小接触压力节点的距离。通过接触压力的分布图,我们找到最大接触压力节点为3141节点,最小接触压力节点为3227节点。我们在命令行中输入b = DISTND(3141,3227),命令输出窗口即显示b的值为0.2488。
对比使用赫兹公式计算出的接触面半宽0.2407mm,ANSYS计算的接触面半宽0.2488mm,误差为3.4%,也是可以接受的。