轮轨蠕滑计算模型的比较研究

轮轨蠕滑计算模型的比较研究

作为车辆动力学的基础理论，f.c.chart首先提出了滑动的概念，并给出了平行同质圆柱之间的纵向滑动力和纵向滑动率之间的非线性关系。1964年K.L.Johnson和P.J.Vermeulen对Carter的工作进行了拓展,提出了接触面是椭圆、蠕滑率任意的纵、横向蠕滑系数计算公式,但仍没有考虑轮踏面的锥度。1967年J.J.Kalker运用三维弹性体理论,使轮轨蠕滑理论提升到了实用阶段,提出了三维(纵向、横向和回旋)蠕滑理论,蠕滑力与蠕滑率的关系是非线性的。其后编制了非线性的DUVOROL计算机程序和简化的FASTSIM计算机程序,以及线性理论的计算公式,在80年代至90年代又推出了采用精确弹性理论的计算机程序CONTACT和USETAB数表。Z.Y.Shen和J.K.Hedrick等人又根据以上理论在1983年提出了实用的计算机公式化近似算法。近年来又提出了复椭圆的非赫兹接触理论和轮缘接触的磨合接触。轮轨蠕滑力与蠕滑率之间的关系是很复杂的,它不仅与蠕滑率有关,还与轮轨之间的接触斑形状、接触压力的分布、轮轨表面状态及轮轨材质等因素有关。因而计算蠕滑力时若要全面考虑这些因素,在工程应用上几乎是不可能的。又由于某些参量的离散性很大,对车辆动力学性能的影响又不敏感,所以精确计算轮轨蠕滑力从工程应用角度考虑也是没有必要的。所以以上各种算法都基于不同的假设,适用于不同的场合。而且这些算法产生于不同的时代,受到了计算技术的影响。随着PC机的迅猛发展,PC机无论在速度上,还是在内存容量上都有了本质的改变。如以前制约DOS程序的数组寻址限制问题在现在的Windows程序上已不再存在,计算速度与最初的PC286相比已提高了几个数量级,内存原来为640KB,现在可以选配128MB以上。因而有必要根据现有的计算条件和工程应用的水平对蠕滑力模型和算法进行分析和评估,以更好地应用于机车车辆的动力学分析计算。1轮轨历史模型模拟两个弹性滚动体在正压力作用下,在相互接触处的切平面上能产生切向力,在该切平面的法线方向能产生回转力矩。这些力和力矩的大小主要与上下接触面之间的蠕滑率等有关,蠕滑率分为切向蠕滑率和回转蠕滑率,切向蠕滑率在两垂直的纵横坐标轴上的投影分别为纵向蠕滑率和横向蠕滑率。切向力在同一坐标轴上的对应分量称为纵向蠕滑力和横向蠕滑力,如图1所示。Fx、Fy和Mz分别为纵向蠕滑力、横向蠕滑力和回旋蠕滑力矩。无量刚化后的各向蠕滑率定义为:vx=Vx′-VxVvy=Vy′-VyVωz=Ωz′-ΩzVvx=Vx′−VxVvy=Vy′−VyVωz=Ωz′−ΩzV蠕滑力除了与这些蠕滑率有关以外还与接触区的形状、正压力的分布、接触区的表面润滑状态、材质特性参数等有关。它们之间的关系一般是非线性的。寻求合理的数学模型,并以工程允许的精度来反映它们之间的关系是轮轨蠕滑模型建模的关键,即在给定蠕滑率和其它参数及边界条件下计算获得具有一定精度的蠕滑力和蠕滑力矩。现在铁路轮轨接触的蠕滑模型主要考虑单点接触,并假定匀质椭球体,即接触区的形状为椭圆。如图1所示,椭圆长轴之半为a,短轴之半为b。由于解决问题的性质不同和发展的渐进过程,常用的模型在形式上可大致分为公式化、复杂程序化和表格化三种。根据模型的细化程度又可以分为线性的、非线性的和准非线性的。1.1旋转活性公式1958年Johnson将Carter的两维理论延伸到两个滚动球体的三维工况,其中包含有纵向蠕滑和横向蠕滑,但没有自旋蠕滑。1964年Johnson与Vermeulen又将光滑的半空间理论引入研究没有自旋蠕滑的纯蠕滑工况。其中提出了对合成蠕滑力的修正公式,但没有计及回旋蠕滑的影响。1967年Kalker提出了计算蠕滑力的线性理论和计算公式,在小蠕滑率下依此计算有一定的精度,且计算方便。对于较大蠕滑率情况下的蠕滑力计算,则结合Kalker的线性理论公式和J-V的修正公式,得到了“缩减因子公式法”(又称沈氏定理),它应归类于公式化的准非线性法。首先用Kalker的线性公式计算线性蠕滑力Fx和Fy:Fx=-G(ab)C11vx,Fy=-G(ab)C22vy-G(ab)3/2C23ωz合成的切向蠕滑力F′r为:F´r=(Fx2+Fy2)12F′r=(Fx2+Fy2)12然后修正合成切向蠕滑力Fr:Fr={μΝ(F´rμΝ)-13(F´rμΝ)2+127(F´rμΝ)3(F´r≤3μΝ时)μΝ(F´r＞3μΝ时)Fr={μN(F′rμN)−13(F′rμN)2+127(F′rμN)3(F′r≤3μN时)μN(F′r＞3μN时)因Kalker的横向蠕滑力计算中包含回转蠕滑的影响,所以在修正公式中也包含了回转蠕滑力矩的成分。但从公式的本身看是不精确的。回转蠕滑力矩Mz的计算公式为:Mz=G(ab)3/2C23vy-G(ab)2C33ωz式中:a、b分别为接触椭圆的长短半轴长;G为刚度模量;μ为摩擦系数;Ν为正压力;C11,C22,C23,C33为蠕滑系数;vx,vy,ωz分别为纵向、横向和自旋蠕滑率。1.2基于高速算法的轮轨求解1978年Kalker根据三维精确理论编制了DUVOROL程序,但算法复杂,计算很慢。1982年编制了CONTACT程序,虽然增强了计算功能,但在计算速度上并没有大的加快。同期Kalker提出基于新简化理论的快速算法程序FASTSIM,为寻求速度,作了较大的简化,但在模型中有错误。Kalker指出由于轮轨表面的状态参数和摩擦系数有很大离散度,对轮轨间的切向力计算的精度要求不必太高。这是在当时的计算机条件下提出的实用途径,因而在80年代FASTSIM程序在曲线通过的计算中得到了广泛应用。1.3计算参数的压缩自适应识别插值表格法最初在1980年被英国Derby铁路研究机构所采用。随着计算机内存容量的大幅度提高和相应软件编译器的发展,Kalker在1995年为了兼顾计算速度和精度又再次推出了插值表格法(USETAB),比英国Derby采用的表格方法精度更高。USETAB由4.5兆的数表TABCON文件和一些辅助程序组成,具体结构如图2所示。TABCON文件内的数据由CONTACT程序计算产生,三列数据分别是作用在轮子上的规格化蠕滑力和蠕滑力矩,是规格化的蠕滑率和a/b的函数。柏松比为0.28。材质的弹性模量G为82000N/mm2。规格化的蠕滑率为:vx′=-G(ab)c11vx3FΝvx′=−G(ab)c11vx3FNvy′=-G(ab)c22vy3FΝvy′=−G(ab)c22vy3FNωz′=-G(ab)32c23ωzFΝωz′=−G(ab)32c23ωzFN规格化的蠕滑力为:Fx′=FxFΝFx′=FxFNFy′=FyFΝFy′=FyFNΜz′=Μz√abFΝMz′=Mzab√FN该方法在1996年由Kalker引入我国,在所有快速算法中是精度最高的。目前在我国尚未推广应用。要将4.5兆的数据读入内存,采用基于DOS的640kB的内存的Fortran编译器是不可能实现的,在微机上可能的途径是用NDPFortran,Fortran90或FortranPowerStation。另一方法是将该程序翻译成其它支持大数组的语言,如VC++等。这里对USETAB作一简单介绍。Kalker罗列了9个与计算蠕滑力有关的变量:vx、vy、ωz、a0/b0、c(√a0*b0)、G、si(柏松比)、μ和正压力N。形成数表时将G、c、Nμ设置为常数1。将柏松比si设置为0.28,得到四个独立的变量:a/b、vx、vy和ωz。将其中三个变量规格化为:vx′=-Gc2C11vx/(3Nμ)vy′=-Gc2C22vy/(3Nμ)ωz′=-Gc3C33ωz/(Nμ)整个表格根据a/b的取值平均分为1/7<a/b≤1和1≤a/b<7两个部分,每个部分又根据s与a/b(在区间上1≤a/b<7上取b/a)的乘积取整后的值平均分成1～7七个单元;再将每个单元按照vx′的取值各分为0<vx′≤1和vx′>1两个部分,然后每个部分再按照s与vx′(vx′>1时取vx′=1/vx′)的乘积取整后的值平均分成0～7八个单元;再将每个单元按照vy′的取值各分为0<vy′≤1和vy′<1两个部分,然后每个部分再按照s与vy′(vy′>1时取vy′=1/vy′)的乘积取整后的值平均分成0～7八个单元;再将每个单元按照ωz′的取值各分为0<ωz′≤1、ωz′>1、0>ωz′>-1、和ωz′≤-1四个部分,然后每个部分再按照s与ωz′(ωz′≥1时取ωz′=1/ωz′;0>ωz′>-1时取ωz′=-ωz′;ωz′≤-1时取ωz′=-1/ωz′)的乘积取整后的值平均分成0～7八个单元。当vx′<0时,有Fx=-Fx(|vx′|),当vy′<0时,也有Fy=-Fy(|vy′|)。因而共有:2×7×2×8×2×8×4×8×3个数据。如表1所示。当a/b大于7或小于1/7时,由于难以用线性插值法,仍用改进后的FASTSIM程序计算。2对于工况参量的影响为便于讨论,以下分别用“J-V”、“FASTSIM”和“USETAB”符号来代表以上三种算法。首先是目标参照对象问题,虽然人们在轮轨蠕滑方面做出了大量的工作,但不同的试验,因条件不一,结果相差很大。理论模型也有多种,在学术届Kalker1982年推出的CONTACT程序被公认为较细化的模型。而USETAB的数据表源于CONTACT程序,它们之间的误差为1.5%。因而这里以USETAB作为相对参照对象。第二个问题是提取信息问题,由于比较的主要数据有纵向蠕滑力、横向蠕滑力和回旋蠕滑力矩,而这些数据与多个参量有关,故是一个多维空间问题。这里将这些参量分为三大类,第一类是主动参量,包括横向蠕滑率和纵向蠕滑率;第二类是工况参量,包括接触面积、形状、摩擦系数、正压力等;第三类是回旋蠕滑率。横向蠕滑率和纵向蠕滑率是可以同时连续变化的参量,比较时必须同时考虑。工况类参量是固定的,比较时可以作为条件参数处理。回旋蠕滑率的大小决定于接触区的锥度和轮径,实际上是可以连续变化的,但对于给定的车轮踏面,由于其外形的锥度变化范围和接触位置几乎是固定的,可以将它处理成三个区间,如在踏面上取较小的值,在轮缘上取较大的值,轮缘根部取中等值。经过这样的处理,横向蠕滑力和纵向蠕滑力的比较问题就成为在不同条件和工况下的两维问题(横向蠕滑率和纵向蠕滑率)。回旋蠕滑力矩的比较问题也可以处理成不同工况下的两维问题(横向蠕滑率和纵向蠕滑率)。关于回旋蠕滑力矩与纵向蠕滑率的关系,由于对称性,可以看成是与合成切向蠕滑率的关系,即可以用单横向蠕滑率来代替。至于取值,可以用与USETAB的绝对误差和相对误差的最大值。关于参量取值的范围问题,这里根据实际可能出现的范围选取。横向蠕滑率和纵向蠕滑率,从零至完全饱和为止,一般为10%;回旋蠕滑率可以根据踏面锥度1/Δ和轮对实际滚动圆半径R求得:ωz=sin(arctg(1/Δ)R1/Δ的变化范围一般为:1/40～2.75(取轮缘角为70°)。对应的Δ为0.36～40,这里取Δ为0.35～40。对于半径R,考虑到不同大小的干线机车车辆车轮和磨耗后的大小变化,R的可能变化范围为:360～625mm。若以360作为实际输入参数,将R的变化合并到Δ的变化当中,则修正后的Δ的变化范围为0.35～70。当计算横向、纵向蠕滑力时,取较小自旋蠕滑率,即Δ取为60。关于工况类参数,这里根据通常出现的情况,选择了以下6种组合,由3种蠕滑水平(全蠕滑、半蠕滑水平和1/4蠕滑水平)和2种接触斑长短轴之比组成。2.1种活性水难以融合小旋流体的活性表2为小自旋蠕滑率情况(即Δ为60)下,3种蠕滑水平时的纵向和横向蠕滑力的差异比较。计算条件是:轮对滚动圆半径为525mm,车轮踏面横断面外形的半径为无穷大,钢轨头部的半径为300mm,接触斑的长短轴之半分别为8.224mm和5.757mm。表3为小自旋蠕滑率情况(即Δ为60)下,三种蠕滑水平时的纵向和横向蠕滑力的差异比较。计算条件是:轮对滚动圆半径为525mm,车轮踏面横断面外形的半径为450mm,钢轨头部的半径为300mm,接触斑的长短轴之半分别为8.224mm和5.757mm。2.2回转园滑率值的计算纵向蠕滑率取为零,对应于不同的横向蠕滑率,当回旋蠕滑率变化时,用USETAB、FASTSIM、J-V三种方法计算所得回旋蠕滑力矩数值的比较结果如表4所示。3轮轨探索的横向动态曲线通过USETAB插值表格法、FASTSIM算法和饱和因子公式算法的比较表明,当轮轨接触点位于踏面上,即回旋蠕滑率较小时,不同算法之间的差异变化比较大。由FASTSIM算法求得