欧拉坐标系下柔性轮对旋转效应对轮轨力的影响

崔 潇，姚建伟，胡晓依，孙丽霞，常崇义

(1.中国铁道科学研究院，北京 100081；2.中国铁道科学研究院集团有限公司 铁道科学技术研究发展中心，北京 100081)

随着世界铁路运输的发展，车辆—轨道耦合动力学得到迅速发展。多体动力学由于其计算速度较快的优势，在车辆—轨道耦合动力学计算中使用广泛。

1996年，Brown[1]对柔性转动轴的惯性力问题进行了研究，基于多体动力学对柔性旋转轴的运动方程进行了重新推导，从而能够计算非恒定角速度的工况。此后Brown和Shabana进一步完成了基于拉格朗日方法的任意截面旋转欧拉梁的计算[2]。Meinders[3]和Meinke[4]基于Brown和Shabana的工作，建立了车辆—轨道动力学模型。将轮对考虑为柔性时可以更真实地反应轮轨实际接触关系，使车辆—轨道耦合动力学研究结果更符合实际[5]。但是这类基于拉格朗日坐标系的柔性轮对模型，由于计算速度的限制，无法考虑柔性轮对的旋转效应。

Baeza[6]在考虑轮对绕轴线转动的前提下，采用了欧拉坐标系表示轮对模态，推导了相应的计算公式，计算了轮对通过扁疤时的车辆—轨道耦合动力学状态。Fayos[7]采用欧拉坐标系表示模态以考虑陀螺效应，并将该方法应用到了圆柱体和轮对上。但均未提出充分考虑柔性轮对旋转效应对其他部件刚性位移计算影响的耦合计算方法。Nielsen[8]使用旋转轮对模型，对接触力/变形量的频率响应函数进行了研究，并对频率响应函数峰值随轮对转动速度变化的规律进行了总结。但仅研究了柔性旋转轮对的固有频率分离特性，尚未综合考虑车辆—轨道整体系统对轮轨力响应频率分离的影响。

本文建立刚柔耦合车辆多体动力学方程，实现欧拉坐标系下旋转效应柔性轮对与拉格朗日坐标系下车辆动力学方程耦合的具体计算方法，计算了欧拉坐标系下的柔性轮对频率响应函数，并对车辆—轨道整体系统下轮轨力响应频率分离现象进行了分析，进而研究了采用旋转效应柔性轮对对车辆—轨道系统轮轨间相互作用力的影响。

1 刚柔耦合车辆—轨道系统多体动力学模型

1.1 拉格朗日坐标系下的多体动力学模型

刚柔耦合车辆—轨道系统多体动力学整体模型主要包括车辆系统和轨道系统，如图1所示。其中，车辆系统包括车体、构架、轮对以及悬挂装置；轨道系统包括钢轨、扣件系统和轨道板等。

图1 刚柔耦合车辆—轨道系统计算模型

模态缩聚后，无约束情况下刚柔耦合车辆—轨道系统多体动力学计算公式[9]为

(1)

式中：下标r表示与刚性相关的量；下标f表示与柔性相关的量；M为质量矩阵；K为刚度矩阵；u为拉格朗日坐标系下的位移矢量；Fv为科氏力与离心力的合力矢量；Fe为广义外力矢量。

对于式(1)，在考虑柔性轮对旋转效应的情况下直接计算，Mfr中柔性部分相关量的计算量较大。因此，将式(1)展开为等价形式式(2)和式(3)，并将式(3)从拉格朗日坐标系转化到欧拉坐标系计算，再结合式(2)求解，可以大幅减少计算量。

(2)

(3)

1.2 欧拉坐标系下旋转效应柔性轮对模型

在欧拉坐标下对旋转效应柔性轮对(既考虑柔性、又考虑旋转)计算公式进行推导时，旋转效应相关量可以采用线性方式进行描述，从而使其在多体动力学在线计算前完成计算，因此能够克服拉格朗日坐标系下计算量大的困难，并能将欧拉坐标系下推导的柔性轮对计算公式，与上一节推导的整体计算公式进行耦合，从而能够较好处理旋转效应柔性轮对的建模问题。

Shabana[2]等将式(3)表示为如式(4)所示的拉格朗日坐标系下的轮对模型，为

(4)

式中：ω为柔性轮对在转动后坐标系中的角速度；W为转角系数矩阵；u0为柔性轮对未变形状态下的位置；θ为柔性轮对转角；N为形函数矩阵；ρ为柔性轮对材料密度；v为车速。

当柔性轮对采用角速度ω=(0Ω0)T时(其中Ω为绕轮对车轴转动的角速度)，式(4)中的Sω，Sωω，Sω0分别简化为

Sω=ΩJ1

(5)

Sωω=-Ω2H1

(6)

Sω0=-Ω2L1

(7)

其中，

式中：J1为拉格朗日坐标系下陀螺效应相关矩阵；H1为拉格朗日坐标系下对流相关矩阵；L1为拉格朗日坐标系下离心力相关矢量；Φ为模态转换矩阵；rL为拉格朗日坐标系下的位置坐标；J0为陀螺效应基本矩阵；H0为对流基本矩阵。

使用正则振型将式(4)中Mff对角化，并将拉格朗日坐标系下的式(4)转换为欧拉坐标系下的式(8)[3]。

(8)

F(rE,t)=α(rE)g(t)

ΦT(rE)=ΦTN(rE)

式中：sE为欧拉坐标系下的模态坐标；G为欧拉坐标系下陀螺效应对偶矩阵；J2为欧拉坐标系下陀螺效应相关矩阵；H2为欧拉坐标系下对流相关矩阵；L2为欧拉坐标系下离心力相关矢量；P为欧拉坐标系下外力相关矢量；Φ(rE)为rE位置处的形函数矩阵；N(rE)为rE位置处的形函数矩阵；F(rE,t)为旋转效应柔性轮对所受的外力，可以表示为2部分的乘积，分别为与位置相关的α(rE)以及与时间相关的g(t)；K为对角阵，且对角线上的元素为轮对固有频率的平方。

J2，H2，P，L2，G无法通过多体动力学商用软件直接计算，采用有限元方法计算得

(9)

(10)

(11)

(12)

(13)

式中：Wt为数值积分权系数；|JF|为雅可比行列式；h为总单元总数；i为单元序号；j为积分点序号。

2 模型求解

以上建立考虑轮对旋转效应的刚柔耦合动力学模型，在求解上存在2个难点：一是如何计算欧拉坐标系下的旋转效应柔性轮对，其计算涉及单元层级，无法通过商用有限元软件直接计算；二是已经解耦的式(2)处于拉格朗日坐标系下，而式(6)处于欧拉坐标系下，如何重新耦合2个处于不同坐标系下的公式，从而进行计算。

2.1 欧拉坐标系下旋转效应柔性轮对求解

采用OpenSees开源有限元软件，添加欧拉坐标系下单元级计算函数，二次开发欧拉坐标下旋转效应柔性轮对的计算功能。主要利用OpenSees软件中计算所涉及的Brick 3维实体单元类[10]和分析类。

在计算过程中，分析类主要起到总体控制的作用，调用FE_Element和Brick类，遍历所有单元，并获取每个单元计算的结果。然后，再对单元计算的结果进行组合，获得全局下的计算结果。FE_Element和Brick类中需要添加2类计算功能：一类主要用于欧拉坐标下计算旋转柔性体所需的J2，H2，P，L2，G矩阵；另一类主要用于计算与柔性体相关的各类矩阵。通过以上的二次开发，能够实现欧拉坐标系下的相关计算。

2.2 欧拉坐标系下旋转效应柔性轮对与整体多体动力学系统耦合求解

通过式(8)计算出的欧拉坐标系下模态坐标sE，耦合求解时需要转换为拉格朗日坐标系下的模态坐标sL，为

sL=BsE

(14)

Bk=

式中：B为欧拉坐标系下柔性坐标与拉格朗日坐标系下模态坐标的转换矩阵，它是1个分块对角矩阵，除了对角线上的分块对角矩阵，其余位置取值都为0；k为模态频率序号；n为选取模态频率的总数；θ为轮对转角。

系统耦合求解的主要流程如下。

(1)获取上一时间步轮对的外力，转化为轮对在欧拉坐标系下受到的弹性力；求解式(8)，计算欧拉坐标系下的模态坐标。

(2)使用式(14)，将欧拉坐标系下的模态坐标转换为拉格朗日坐标系下的模态坐标。

(3)将转换后的模态坐标结合式(2)进行求解。计算轮对所受的外力，进入下一时间步的流程(1)。

通过将轮对所受外力转化到欧拉坐标系，将旋转轮对的欧拉坐标系下的模态坐标转化到拉格朗日坐下的模态坐标，不仅可以使欧拉坐标系下的旋转效应柔性轮对计算与拉格朗日坐标系下的整体多体动力学系统计算实现耦合，而且由于每一次计算都使用的是当前步的模态坐标，能够在使用式(2)计算时，充分考虑轮对的柔性、旋转效应，从而实现迭代计算。

3 计算分析

3.1 车辆—轨道系统参数

针对我国高速铁路车辆—轨道系统，车辆模型参数采用我国某高速动车组基本参数，其中车轮直径为0.86 m，踏面为LMA型；轨道模型参数选用某类型无砟轨道参数，轨枕间距为0.65 m，钢轨廓形为CHN60，设1/40轨底坡。车辆—轨道系统模型部分参数见表1。为了避免钢轨两端引起的应力波效应，且考虑轨道扣件整体刚度及阻尼作用[11]，建立的轨道全长为100 m，车辆从40 m处开始运行，设运行初始位置为里程零点，共运行里程约为12 m。

表1 车辆—轨道系统模型部分参数

3.2 轮轨力频率特性

选取车速为300 km·h-1，分别建立刚性轮对、无旋转效应柔性轮对以及旋转效应柔性轮对模型。柔性轮对的旋转效应是由柔性和旋转共同引起的。刚性轮对只考虑轮对的旋转，但并不考虑其柔性；无旋转效应柔性轮对只考虑了柔性，但不考虑轮对绕车轴轴线的旋转；旋转效应柔性轮对既考虑了柔性又考虑了旋转。对于刚性轮对和无旋转效应柔性轮对，采用拉格朗日坐标系下的模型；旋转效应柔性轮对则采用欧拉坐标系下的模型。

对比分析这3种轮对建模条件下列车通过直线轨道时轮轨间相互作用情况。1位轮对轮轨间相互作用力对比结果如图2所示。

图2 不同轮对建模条件下轮轨作用力

从图2可以看出：在不同的轮对建模条件下，轮轨间作用力存在明显差异；由于扣件系统的离散支撑作用，轨道存在固有的刚度不平顺，轮轨间作用力随之出现波动，刚性轮对通过直线轨道时，轮轨间作用力波动较小，主要存在单一频率的周期性波动，其中轮轨垂向力平均值约为60 kN，波动幅值约为0.05 kN；无旋转效应柔性轮对条件下，轮轨垂向力出现等幅周期性波动，幅值约为1 kN，并存在1个明显的高频周期性波动，幅值约0.3 kN；旋转效应柔性轮对条件下，轮轨垂向力波动情况相对复杂，存在多个主频特征，且比无旋转效应柔性轮对的高频波动更加剧烈；轮轨横向力波动情况与轮轨垂向力类似；刚性轮对的轮轨纵向力基本为0；考虑旋转效应的柔性轮对波动明显大于不考虑旋转效应时，两者存在不同相位，且无旋转效应柔性轮对的波动幅值基本为定值，而旋转效应柔性轮对的波动幅值没有明显规律。

为了进一步研究轮轨作用力的波动规律，以轮轨垂向力为例，对这3种轮对建模条件下轮轨垂向力做频谱分析，输出频率为5000 Hz。取车辆稳定运行后里程为4～12 m的轮轨垂向力进行频域分析，如图3所示，频率范围取20～2 000 Hz。

图3 轮轨垂向力频域分析

柔性轮对本身的固有模态也对轮轨力有较大影响，因此需对无旋转效应柔性轮对进行模态分析。其1，3，5阶弯曲模态频率分别为124.3，393.5，863.3 Hz，其振型如图4所示。

从图3可以看出如下结果。

(1)刚性轮对的轮轨垂向力存在约128.2 Hz的主频，通过计算可知，该频率与列车以300 km·h-1速度通过间距为0.65 m轨枕的激励频率相一致，且功率谱密度较小，该主频不仅存在于刚性轮对的轮轨垂向力中，也存在于柔性轮对的轮轨垂向力中，但由于所占成分较少，在柔性轮对的功率谱中并不明显。

图4 无旋转效应柔性轮对1，3，5阶弯曲模态

(2)无旋转效应柔性轮对的轮轨垂向力存在3个主频，分别为121.1，391.7，861.8 Hz，其中121.1 Hz主频对应的功率谱密度较大，这3个主频分别与无旋转效应柔性轮对的3个弯曲模态频率接近(如图4所示)。由此可知，无旋转效应柔性轮对的轮轨垂向力主要出现了这3个模态频率的叠加振动，即无旋转效应柔性轮对通过直线轨道时，将不同程度地激起轮对的1，3，5阶弯曲振动。

(3)旋转效应柔性轮对的轮轨垂向力存在多个主频，其中35.61 Hz低频接近列车以300 km·h-1速度运行时轮对转动1周的频率，92.59和149.6 Hz频率接近轮对1阶弯曲模态频率，377.5和420.2 Hz频率接近轮对3阶弯曲模态频率，833.3和890.3 Hz频率接近轮对5阶弯曲模态频率。

对于旋转效应柔性轮对而言，其固有模态特性还与其旋转速度相关，可通过频率响应函数进行描述。分析时分别选取柔性轮对不旋转以及旋转速度为193.8 rad·s-1即车速为300 km·h-1这2种工况。在车轮踏面名义滚动圆位置输入垂向激励，频率范围为0～2 500 Hz，获取轮辋上某点垂向位移的响应，输出频率为5 000 Hz，其0～1 200 Hz频段的频率响应函数如图5所示。

根据图5可以看出：速度为0 km·h-1时，旋转效应柔性轮对的弯曲模态频率与无旋转效应柔性轮对的一致，随着速度的增加，旋转效应柔性轮对的弯曲模态频率由1个分离为2个。由此推断，由于旋转效应柔性轮对发生了频率分离现象，轮轨垂向力的频率也发生了频率分离，且分离后频率1个略大于分离前频率、1个略小于分离前频率；轮轨垂向力响应的1 040 Hz频率接近轮对1 043 Hz的固有频率，且功率谱密度较大，表明旋转效应柔性轮对与钢轨某固有频率发生了共振。

图5 柔性轮对频率响应函数

为进一步研究车速对旋转效应柔性轮对轮轨垂向力频率分离现象的影响，利用刚柔耦合车辆—轨道系统动力学模型分析刚体轮对、无旋转效应柔性轮对和旋转效应柔性轮对的主频特征。不同车速下3种轮对的轮轨垂向力主频如图6所示，图中f1，f3和f5分别为无旋转效应柔性轮对1，3，5阶弯曲模态频率；f11和f12分别为旋转效应旋转效应柔性轮对1阶弯曲模态频率分离后的2个频率，f31和f32，f51和f52分别对应3，5阶弯曲模态频率分离后的2个频率。

图6 轮轨垂向力主频

从图6可以看出：刚性轮对的轮轨垂向力主频与轨枕间距直接相关，当轨枕间距为定值时，其随车速线性增加；无旋转效应柔性轮对的轮轨垂向力主频不随车速改变；旋转效应柔性轮对的轮轨垂向力随车速增加时，其主频分离现象更加明显，分离后这2个频率的差值随车速的增加而增大；当车速大于350 km·h-1时，1，3，5阶弯曲模态频率的分离频率差值在50 Hz以上。

因此，在研究车辆—轨道系统高速轮轨动态响应问题时，采用考虑旋转特性的柔性轮对格外重要。此外，本节建立的轨道模型并未存在轨道几何不平顺，但当旋转效应柔性轮对通过该区段时仍引起了多个固有频率振动的叠加，表明对于新建或新修高速线路，柔性轮对旋转效应易引起直线轨道与轮对1，3，5阶弯曲模态分离后频率相对应的轮轨垂向力。

3.3 周期性轨道不平顺激励下轮轨力响应

由前文可知，不同于刚性轮对和无旋转效应柔性轮对，旋转效应柔性轮对轮轨垂向力在1 040 Hz频率附近出现了功率谱密度较大的主频，该主频与轮对1 043 Hz的固有频率接近。为了进一步研究该频率对车辆—轨道系统高速轮轨力的影响，建立与该频率相对应的周期性轨道不平顺模型，以车速300 km·h-1为例，轨道不平顺波长设为80 mm。考虑到钢轨表面出现0.02 mm周期性轨道不平顺时需要进行打磨，因此将轨道不平顺设为0.02 mm幅值的简谐波，从里程5 m处开始设置。同时，考虑在轨道平顺区段轮轨力存在一定波动(见图2)，因此，将轨道不平顺和平顺区段的结果进行对比，获得3种轮对建模条件下轨道不平顺区段轮轨作用力，如图7所示。

为了更加清楚地说明柔性旋转效应的影响，计算周期性轨道不平顺条件下无旋转效应柔性轮对和旋转效应柔性轮对的轮轨力与刚体轮对之间的差值，分别用“柔—刚”和“柔转—刚”表示，结果如图8所示。

结合图7和图8可以看出：在轨道不平顺区段，轮轨力随轨枕位置出现周期性波动，在轨枕附近轮轨力振动较大，在2轨枕间振动较小，轮轨力振动幅值与轨枕位置之间存在相位差，其中轮轨垂向力波动对轨枕附近钢轨冲击较大，最大幅值约为40 kN；柔性轮对旋转效应对轮轨垂向力和横向力的影响不明显，考虑旋转效应时与刚体轮对的差值略大于忽略旋转效应；轮轨纵向力受柔性轮对旋转效应影响较明显，特别是当轨道出现短波不平顺时，旋转效应柔性轮对的轮轨纵向力高频振动更剧烈，在波长为80 mm的轨道不平顺作用下，相比于刚性轮对基本为零的轮轨纵向力和无旋转效应柔性轮对轮轨纵向力约0.1 kN的振幅，旋转效应柔性轮对轮轨纵向力的最大波动峰峰值到达了1.8 kN。因此，对于研究高速铁路系统轮轨短波病害问题时，轮对建模应考虑其旋转效应。

图7 周期性轨道不平顺区段轮轨作用力

为了研究柔性轮对轮轨动态响应对应的主频对车辆—轨道系统轮轨力的影响，分别建立与轮对1，3，5阶弯曲模态频率和1 043 Hz主频相对应的周期性轨道不平顺模型，其波长分别为670.0，211.8，96.5和80.0 mm；轮对模型既考虑柔性，又考虑旋转，车速仍设为300 km·h-1，得到不同主频下旋转效应柔性轮对轮轨力，如图9所示。

由图9可以看出：随着轨道几何不平顺波长的减小，列车通过不平顺区段频率的增加，轮轨垂向力和横向力的波动呈上升趋势，而轮轨纵向力波动的趋势不明显，其主要随轨枕位置波动。值得注意的是：在80 mm波长轨道几何不平顺情况下，轮轨垂向力和横向力出现了随轨枕位置里程波动的时大时小的“拍振”现象，在两轨枕间轮轨力的波动甚至小于96.5 mm时的情况。结合图3(c)无不平顺区段的频域分析可知，此“拍振”现象是由于车辆—轨道系统出现了共振。对于高速铁路线路应避免该频率对应的轨道短波周期性不平顺的出现(如波磨等)，如出现应尽早进行轨道相关养护维修作业。

图8 周期性轨道不平顺和平顺区段轮轨力差值对比

4 结 论

(1)通过将轮对所受外力转化到欧拉坐标系，并将旋转轮对在欧拉坐标系下的模态坐标转化到拉格朗日坐标系下，可以实现欧拉坐标系下的旋转效应柔性轮对计算与拉格朗日坐标系下的车辆—轨道多体动力学整体模型计算之间的耦合，进而完成迭代计算。

图9 不同主频下旋转效应柔性轮对的轮轨力

(2)考虑旋转效应时，柔性轮对固有频率出现了频率分离现象，且随着轮对旋转速度的增加，频率分离现象越加明显，分离频率大小值之间差距越大。当车速大于350 km·h-1时，与旋转效应柔性轮对的1，3，5阶弯曲模态相对应的轮轨垂向力分离频率差值将大于50 Hz。

(3)在扣件系统离散支撑作用下，轨道存在刚度不平顺，轮轨间作用力随之出现波动。柔性轮对的轮轨间垂向力基本为等幅周期性波动，而旋转效应柔性轮对轮轨间垂向力波动情况相对复杂，存在多个主频特征，且比无旋转效应柔性轮对的高频波动更加剧烈。

(4)旋转效应对柔性轮对轮轨纵向力的影响较垂向力和横向力明显，当高铁线路存在轨道几何短波不平顺时，旋转效应易加剧柔性轮对纵向力高频波动；对于轨枕间距为0.65 m的高铁线路，当列车在波磨区段的通过频率约为1 043 Hz时，轮轨间垂向力出现了随轨枕位置波动的“拍振”现象，即引起了车辆—轨道系统共振，应尽早对其进行轨道养护维修作业。