一种简化轮轨非赫兹切向接触计算的等效方法

1.本发明涉及铁道工程轮轨接触计算领域，具体是一种简化轮轨非赫兹切向接触计算的等效方法。


背景技术：

2.轮轨系统是连接车辆与轨道的重要纽带，其中，宏观的轮轨蠕滑力是计算车辆-轨道耦合动力学的重要参量，轮轨滚动接触力学行为则为求解轮轨滚动接触疲劳和磨损问题提供必要的基础解。目前在铁路行业中得到广泛应用的接触模型主要有fastsim，沈氏理论与polach方法以及contact。然而，这些模型仍存在两点不足：一方面，如fastsim，沈氏理论和polach方法等简化模型均是基于赫兹假设建立，无法更为准确地描述轮轨非赫兹接触斑；另一方面，如contact等准确的计算模型需要大量的迭代过程，其计算效率难以满足如车辆动力学、轮轨磨损预测等大规模计算要求。
3.现有技术中较高效的模型为stripes模型，stripes模型采用的是局部椭圆法。stripes使用改进的fastsim来进行轮轨切向接触分析。首先将接触斑沿纵向划分为若干条带，并假设每个条带位于局部椭圆的中心，如图1所示，stripes模型针对每个条带计算出局部椭圆，基于局部椭圆计算每个条带上局部的kalker线性蠕滑系数cij和fastsim柔度系数，其求解局部椭圆的步骤如下：
4.首先通过法向接触得到每个条带对应的局部接触椭圆长半轴a(y):式中，g(y)为轮轨间隙函数，a(y)为修正后的轮轨纵向组合曲率。根据轮轨实际型面求解轮轨横向组合曲率式中，cw和cr分别是车轮和钢轨型面的横向曲率。
5.b(y)需要进行平滑和修正，平滑是由于某些位置的轮轨型面曲率半径变化剧烈，修正的原因是根据轮轨型面曲率得到的b(y)有可能为负值，但是根据hertz方法中，b(y)必须恒为正，因此需要将b(y)中的负值替换掉。根据a(y)与b(y)求解局部接触椭圆短半轴：式中δ为刚性渗透量，ε为渗透量修正系数，依据刚性接触点处的轮轨曲率确定：式中，n和r由赫兹理论计算得到。依据每个条带的a与b值，可得到每个条带对应的kalker线性蠕滑系数cij和fastsim柔度系数，进而求解非赫兹工况下的切应力。
6.在stripes模型中，求解a(y)之前，需要对潜在接触区域的轮轨横向组合曲率b(y)进行平滑与修正。主要原因如下：首先，根据hertz理论，b(y)应该恒为正值。对于大多数轮轨接触工况，通过轮轨型面实际曲率计算得到的b(y)为负值，此时需要把负曲率变为正值，否则不能应用于基于hertz理论的切向接触模型。其次，一般来讲，接触斑内轮轨型面曲率半径变化剧烈，而接触斑形状在曲率半径不连续点处会出现突变。实际上，接触斑形状应该
是光滑的，因此需要对其横向曲率进行人为的平滑处理。而平滑方式的选择较主观，且接触结果易受平滑方式的影响，有文献采用移动平滑方法对磨耗的lma车轮踏面曲率进行平滑，对比了不同平滑窗口的移动平滑方法对轮轨法向力的影响，如图2所示。由此可以发现平滑方式的选择会明显改变轮轨接触力。在工程计算时，很难评判哪种平滑方式最优，导致stripes模型稳定性不足，计算精度依赖于人为的主观干预。因此，以往研究中通过修正轮轨型面曲率的方法的计算稳定性有待提高。


技术实现要素：

7.本发明目的是针对背景技术中存在的问题，提出一种将非椭圆接触斑等效为若干局部椭圆的方法，可将基于hertz假设建立的切向接触模型应用于非赫兹接触工况，从而大幅度提升非赫兹切向接触的计算效率的简化轮轨非赫兹切向接触计算的等效方法。
8.本发明的技术方案：一种简化轮轨非赫兹切向接触计算的等效方法，包括以下步骤：
9.s1、轮轨接触斑压应力分布p0(y)、纵向接触边界a(y)和车轮在接触点处的滚动半径r作为已知条件给定，假设轮轨接触满足弹性半空间假设，则纵向曲率在所有条带保持恒定，将接触斑划分为若干沿滚动方向的条带，假设每个接触条带位于局部椭圆的中心，为了在每一个条带上应用切向接触算法，需要求解接触斑内横向曲率分布b(y)和每个局部椭圆的横向半轴长度b(y)；
10.s2、根据hertz理论，如果a＞b，公式一，公式二，如果a＜b，公式三，公式四,公式四,上述公式中，k(e)和e(e)分别是第一类椭圆积分和第二类椭圆积分；获得e2后，椭圆接触斑横向半轴长度b和横向组合曲率b可进行求解；但不能提前知道a与b的大小关系，发明trial算法求解b和b；
11.s3、定义两个函数，公式五，公式六，
12.s4、定义一个极小值c＝1e-9
，得到h1与h2，公式七，h1＝f1(c)
·
f1(1-c)，公式八，h2＝f2(c)
·
f2(1-c)；
13.s5、如果h1＜0，则a＞b，e的值可通过f1(e)＝0得到，b由公式二求得；如果h2＜0，则a＜b，e的值可通过f2(e)＝0得到，b由公式四求得；
14.s6、此外，存在一种特殊的情况e＝1，b＝0，即共形接触，由于无法通过椭圆积分求解该种接触类型，因而其不包含于s5，因此，设定b＝10a。
15.与现有技术相比，本发明具有如下有益的技术效果：提出一种局部椭圆算法trial，利用法向应力分布和接触斑边界计算得到每一横向位置处的横向曲率与每个条带的等效椭圆宽度b，可将基于赫兹假设的切向接触算法(如fastsim)应用于非赫兹工况，从而大幅度提升非赫兹切向接触的计算效率。且利用所提出的方法计算求得局部椭圆，能够避免由轮轨型面获得的横向组合曲率b(y)出现负值或不平滑的情况，即克服stripes方法的局限性。
附图说明
16.图1为stripes模型求解等效椭圆示意图；
17.图2为横向曲率平滑对法向力的影响的示意图；
18.图3为trial方法流程图；
19.图4为trial与通过轮轨型面计算得到的横向组合曲率对比示意图；
20.图5为旋转体示意图；
21.图6为双半椭圆接触的法向解示意图；
22.图7为切向接触解示意图。
具体实施方式
23.实施例一
24.如图1-7所示，本发明提出的一种简化轮轨非赫兹切向接触计算的等效方法，包括以下步骤：
25.s1、轮轨接触斑压应力分布p0(y)、纵向接触边界a(y)和车轮在接触点处的滚动半径r作为已知条件给定，假设轮轨接触满足弹性半空间假设，则纵向曲率在所有条带保持恒定，将接触斑划分为若干沿滚动方向的条带，假设每个接触条带位于局部椭圆的中心，为了在每一个条带上应用切向接触算法，需要求解接触斑内横向曲率分布b(y)和横向半轴长度b(y)；
26.s2、假设两个接触体的材料特性相等，其弹性模量和泊松比分别为e和v，轮轨椭圆接触斑长半轴为a，最大应力为p0；根据hertz理论，，求解局部接触椭圆的短半轴b或横向曲率b；如果a＞b，公式一，公式二，公式二，如果a＜b，如果a＜b，公式四,公式四,上述公式中，k(e)和e(e)分别是第一类椭圆积分和第二类椭圆积分；获得e2后，椭圆接触斑横向半轴长度b和横向组合曲率b可进行求解；但不能提前知道a与b的大小关系，发明trial算法求解b和b；
27.s3、定义两个函数，公式五，公式六，
28.s4、定义一个极小值c＝1e-9
，得到h1与h2，公式七，h1＝f1(c)
·
f1(1-c)，公式八，h2＝f2(c)
·
f2(1-c)；
29.s5、如果h1＜0，则a＞b，e的值可通过f1(e)＝0得到，b由公式二求得；如果h2＜0，则a＜b，e的值可通过f2(e)＝0得到，b由公式四求得；
30.s6、此外，存在一种特殊的情况e＝1，b＝0，即共形接触，由于无法通过椭圆积分求解该种接触类型，因而其不包含于s5，因此，设定b＝10a。
31.本实施例提出一种局部椭圆算法trial，利用法向应力分布和接触斑边界计算得到每一横向位置处的横向曲率，可将基于赫兹假设的切向接触算法(如fastsim)应用于非赫兹工况，从而大幅度提升非赫兹切向接触的计算效率。且利用所提出的方法计算求得局部椭圆，能够避免由轮轨型面获得的横向组合曲率b(y)出现负值或不平滑的情况，即克服stripes方法的局限性。
32.以lmb车轮踏面与chn60钢轨型面为例，钢轨轨底坡设为1：40,车轮半径为460mm。图4对比了通过本方法计算得到的横向组合曲率b以及由轮轨型面得到的结果，其中轮对横移量为0与6mm时。可以发现，1、由轮轨型面得到的横向组合曲率出现了负值，导致该结果无法与基于hertz假设的切向接触模型联合使用；2、在接触点附近，由轮轨型面获得的曲率结果存在突变，这会导致接触结果不准确。相较而言，本发明提出的trial方法改善了此问题，得到的横向组合曲率变化较为平缓且数值未出现负值，不需要再进行修正即可直接应用到接触模型。
33.实施例二
34.如图1-7所示，列举一个利用本发明提出的一种简化轮轨非赫兹切向接触计算的等效方法求解局部接触椭圆进而得到切向接触解的应用案例，具体步骤如下：
35.s1、假设一个旋转体与平面发生接触，该旋转体的横向曲率半径在接触点附近由100mm变化为500mm，旋转体的滚动圆半径为460mm。以kalker的contact程序所得结果作为参考。定义接触坐标系，如图5所示，原点o为接触点，x轴为纵向，z轴垂直于接触平面向下，y轴通过右手法则确定，垂直于x轴与z轴，指向右侧。
36.s2、已知轮轨接触点(0,0)和滚动圆半径r＝460mm，通过contact法向算法得到轮轨接触斑每个条带对应的接触斑纵向边界a(y)，以及法向应力分布p0(y)，结果如图6(a)和(b)所示，纵向曲率a＝1/2/r，r为接触点处的滚动圆半径。
37.s3、采用trial方法，根据每个接触条带计算局部接触椭圆，首先定义两个函数，s3、采用trial方法，根据每个接触条带计算局部接触椭圆，首先定义两个函数，
38.s4、定义c＝1e-9
，得到h1与h2，h1＝f1(c)
·
f1(1-c)，h2＝f2(c)
·
f2(1-c)。
39.s5、如果h1＜0，则e2数值范围为[c,1-c]，则a＞b，e的值可通过f1(e)＝
0得到，求得e之后，通过下式求解接触斑横向半轴长度b(y)与横向组合曲率b(y)，如果h2＜0，则e2数值范围为[c,1-c]，则a＜b，e的值可通过f2(e)＝0得到，求得e之后，通过下式求解接触斑横向半轴长度b(y)与横向组合曲率b(y)，e＝1时，不能采用椭圆积分，它存在于共形接触工况，此时b＝0，这种情况不能采用基于hertz的切向接触模型，因此，引入一个假设，b(y)＝10a(y)，b(y)＝c。
[0040]
s6、基于局部椭圆计算每个条带上局部的kalker线性蠕滑系数cij，cij是关于接触斑长短半轴a,b的比值的函数。
[0041]
s7、每个接触斑条带所对应的局部柔度系数为，s7、每个接触斑条带所对应的局部柔度系数为，其中k(y)＝a(y)/a0。
[0042]
s8、在每个条带上应用fastsim模型，结果如图7所示，可以发现通过trial fastsim得到的切应力与contact结果吻合。
[0043]
上面结合附图对本发明的实施方式作了详细说明，但是本发明并不限于此，在所属技术领域的技术人员所具备的知识范围内，在不脱离本发明宗旨的前提下还可以作出各种变化。