基于ANSYS二次开发的轮轨耦合模型及应用研究

孟宪金，刘林芽，秦佳良，李 辉，左志远

（华东交通大学 铁路环境振动与噪声教育部工程研究中心，南昌 330013）

车辆-轨道相互作用一直是车辆动力学和轨道结构振动研究中最重要的基础性问题[1]。近年来，随着我国高速铁路的快速发展，人们越来越关注这个问题。早期研究并未实现车辆系统和轨道系统的耦合，并且车辆系统和轨道系统的动力响应采用分离的方法计算[2－4]，其结果误差较大，不能满足实际要求。事实上，车辆运行导致轨道振动，这反过来加剧了车辆振动。因此，通过轮轨相互作用关系，车辆系统和轨道系统本质上形成了一个完整的动力学系统。近年来，许多学者以此为基础建立了多种车辆-轨道耦合模型，这些模型在铁路工程中也被广泛采用[5－8]。翟婉明等[9]建立的车辆-轨道耦合动力学统一模型，将显式积分法和模态叠加法应用到无砟轨道动力学研究中，能够实现车辆和轨道运动方程的耦合，计算精度和稳定性都比较高。张斌等[10]借助有限元法和交叉迭代数值算法建立了车辆-轨道-路基非线性耦合系统振动模型，此方法在加速系统迭代的同时，增强了迭代的稳定性，能够较好反映轮轨相互作用引起的系统动力学响应特征，但由于在每一个积分步长内一般需要迭代计算多次，相较于上述方法计算效率略低。在模型建立方面，朱志辉等[11]运用虚拟激励法和有限元法相结合的方法，建立了车辆-轨道-桥梁垂向耦合系统随机振动模型，分析了简支箱梁桥墩顶在轨道不平顺种类和车速影响下的垂向动反力随机特性，但模型中的轮轨接触关系采用Hertz线性接触，这可能导致非线性结构和系统的计算产生误差；闫斌等[12]为探讨桥上无砟轨道损伤对列车-轨道-桥梁系统动力响应的影响规律，基于车辆-轨道-桥梁耦合动力学原理，通过ANSYS 和SIMPACK 联合仿真，建立了双线32m简支箱梁桥CRTSⅢ型无砟轨道空间动力学模型,此建模方法虽能够精确仿真车辆模型，但建模及计算方法严重依赖系统构型。如今，车辆-轨道相互作用理论已经成熟和完善，如何使轮轨耦合相互作用的计算更加方便且更具适用性已成为研究中的关键问题。

本文在以往研究的基础上，基于ANSYS的二次开发技术，利用APDL（ANSYS Parameter Design Language）语言提出一种具有广泛适用性的轮轨耦合相互作用计算方法，并运用此计算方法建立车辆-轨道垂向耦合模型。其中，车辆部分采用多体动力学理论建模，轨道部分采用有限元理论仿真，车辆系统和轨道系统分别采用显式积分和隐式积分求解。然后，通过与文献[13]中的车辆-轨道耦合系统交叉迭代算法的计算结果进行对比，来验证模型和计算方法的正确性。在此基础上，以高速列车-CRTSⅡ型板式无砟轨道为例，利用该方法对不同扣件失效数量工况下车辆-轨道垂向耦合系统的动力响应进行研究，为轮轨相互作用的研究和轨道结构的设计提供理论基础。

1 车辆-轨道垂向耦合振动模型

APDL 作为ANSYS 中的参数化设计语言，拓宽了ANSYS的功能和应用范围，便于建模和仿真。采用这种语言，可以在ANSYS计算平台上轻松考虑车辆系统和轨道不平顺。因此，本文基于ANSYS二次开发语言APDL研究轮轨耦合相互作用。以高速列车-板式无砟轨道为例，建立车辆-轨道垂向耦合模型（如图1所示），其中车辆系统采用多体动力学理论建模，轨道系统采用有限元理论建模，通过轮轨相互作用关系实现车辆系统和轨道系统的耦合。

图1 车辆-轨道垂向耦合模型

1.1 车辆模型

列车系统中，每辆车被模拟成一个以速度v运行于轨道结构上的多刚体半车模型，包括一个车体，两个转向架构架，四个轮对和两级悬挂。在车辆模型中，考虑了车体的沉浮（Zc）运动和点头（θc）运动，前后转向架构架的沉浮（Zt）运动和点头（θt）运动，以及4 个轮对的垂向振动（Zwi，i=1～4），共10 个自由度。具体的车辆动力学方程见参考文献[12]。此外，列车系统的动力学方程可以用矩阵形式表示为：

式中：M、C和K分别为车辆的质量、阻尼和刚度矩阵；Z为车辆的位移向量，包括车体垂向位移和点头位移、转向架的垂向位移和点头位移以及4 个轮对的垂向位移；F为轮轨接触力向量。

此外，式(1)可简化为：

式中：M-1为车辆质量矩阵的逆矩阵。

对于车辆的质量、刚度和阻尼矩阵，无法通过ANSYS 自动计算，应通过APDL 语言手动写入ANSYS中。

1.2 轨道模型

在车辆-轨道垂向耦合振动模型中，轨道由钢轨、轨道板、底座板和路基组成。假定整个系统模型沿线路方向左右对称，因此，轨道子系统取轨道对称轴的一半建模。轨道系统中，采用Beam188 梁单元模拟钢轨、轨道板和底座板，采用Combin14 线性弹簧单元模拟扣件系统和CA砂浆层。

1.3 轮轨相互作用关系

车辆-轨道耦合系统模型属于非线性问题研究，而耦合系统中作为连接车辆子系统和轨道子系统纽带的轮轨相互作用关系，也应是非线性的。因此，轮轨相互作用关系可用非线性Hertz接触理论描述，从而确定轮轨之间的垂向作用力：

式中：G为轮轨接触常数（m/N2/3）；δZ(t)为轮轨间的弹性压缩量（m）。

对于锥型踏面车轮：

对于磨耗型踏面车轮：

式中：R为车轮半径（m）。

在式(3)中，轮轨间的弹性压缩量δZ(t)的表达式为：

式中：Zw(t)为t时刻车轮的位移；Zr(t)为t时刻车轮位置处的钢轨位移；Z0(t)为t时刻车轮位置处钢轨的位移不平顺。

由于钢轨采用有限元法建模，在每个计算步骤中，车轮的位置在多数情况下并不处在钢轨节点位置，而是位于两个相邻节点之间。因此，有必要将轮轨力分配到相邻两节点上，如图2所示。将轮轨力P0施加到梁单元上的第i和第i+1 个节点之间的点上。

图2 钢轨梁单元上轮轨力的分配

进一步，梁单元上的轮轨力可以分配为相邻两节点的两个剪力和两个弯矩，其分配公式可表示为：

式中：Pi和Pi+1分别为梁单元上第i和第i+1 个节点上分配的剪力；Mi和Mi+1分别为第i和第i+1个节点上分配的弯矩；l为梁单元长度；a和b分别为车轮位置距梁单元上第i和第i+1个节点的距离。将轮轨力按式(7)分配后，即可将轮轨力施加到钢轨有限元模型上。

此外，由于钢轨有限元模型中只得到了各节点处的振动响应，为了得到任意相邻两节点之间车轮位置处的钢轨位移，需对相邻节点的振动响应进行处理。本文计算了模型中所有节点的位移，并基于插值理论对任意相邻两节点的振动响应进行整合，得到轮轨力作用点处的钢轨位移，其整合公式可表示为：

式中：Zr为车轮位置处的钢轨位移；Zi和Zi+1分别为第i和第i+1个节点的位移；Ri和Ri+1分别为第i和第i+1个节点的转角；S1、S2、S3和S4为系数，可以用Hermite形状函数表示为：

因此，结合式（7）、式（8）和式（9）可以方便地将轮轨力分配到相邻的两节点上，得到任意位置的钢轨位移，从而确定轨道结构的动力响应。

1.4 数值积分方法

在本文建立的车辆-轨道垂向耦合模型中，车辆系统采用多体动力学理论建模，轨道及下部结构采用有限元法建模，然后利用APDL二次开发语言，将车辆系统方程和轮轨相互作用关系编程到ANSYS中。因此，在数值计算中，采用显式积分和隐式积分相结合的方法对该动力学系统进行计算。

对于车辆模型以及轮轨相互作用关系，为了提高其计算效率，采用一种显式积分法（翟方法）求解。翟方法属于显式二步数值积分法，其基本原理是：下一步的位移和加速度值通过前一步和当前步的位移、速度和加速度值确定，然后根据所列的系统方程求解得到下一步的加速度值，如此循环递推，其积分格式为[14]：

式中：∆t为时间积分步长；下标n-1、n和n+1分别代表第n-1、n和n+1 个子步；ψ、φ为积分常数，研究表明：当ψ=φ=0.5时，显式积分法具有良好的稳定性和较高的精度，因此这里的ψ和φ均取0.5。

对于轨道有限元模型，采用Newmark-β隐式积分法求解，其格式为[15]：

式中：∆t为时间积分步长；下标n、n+1 分别代表第n和n+1 个子步；β、γ是Newmark-β积分控制参数，其值由积分稳定要求确定，基于此积分方法，便可计算得到钢轨及下部结构有限元模型的振动响应。

此外，由于采用两种积分方法计算，为了保证计算的收敛性以及计算数据的稳定性，需选取相同积分步长进行计算，本数值计算中积分步长取1×10-5s。

2 模型验证

为了研究高速行车下车辆-轨道垂向耦合系统时域振动响应，本文采用中国高速铁路无砟轨道不平顺谱，其高低不平顺谱密度的一个时域样本，如图3所示。

图3 中国高速铁路无砟轨道高低不平顺谱

为了验证模型和计算方法的正确性，计算高速铁路列车以300 km/h的速度通过板式无砟轨道时车辆和轨道结构的动力响应，并与文献[13]中的车辆-轨道耦合系统交叉迭代算法的计算结果对比。车辆参数见表1，轨道参数见表2。车辆-轨道系统动力响应计算结果如图4所示。

表1 CRH380高速客车动力学参数

表2 轨道结构参数

从图4中车辆-轨道系统动力响应对比结果可以看出，总体上，本文计算得到的钢轨垂向位移与文献[13]中采用交叉迭代算法得到的结果基本吻合，证明本文的车辆-轨道垂向耦合模型和计算方法是正确的。

图4 车辆-轨道系统动力响应对比

3 算例分析

本文利用该方法分析轨道系统中扣件局部失效情况下车辆-轨道垂向耦合系统的动力响应。数值算例中，车辆系统采用CRH380型动车组，轨道系统采用CRTSⅡ型板式无砟轨道，轨道不平顺采用中国高速铁路无砟轨道不平顺谱，模拟轨道总长度为150 m，扣件类型为WJ-7型，扣件间距为0.6 m，共计251 个扣件。针对轨道系统中可能出现的扣件局部失效情况，扣件按位置顺序进行编号，选取钢轨跨中位置扣件及其左右相邻两个扣件作为研究对象，分析扣件不同失效数量对系统各结构动力响应的影响。

3.1 单个扣件失效对系统振动特性的影响

首先考虑单个扣件失效情况，假设钢轨跨中位置的扣件失效，此时失效扣件的刚度和阻尼系数均为0。为避免其他因素对计算结果的影响，计算时考虑的轨道不平顺为同一不平顺。假设列车的运行速度为300 km/h，分别计算扣件失效前后车轨系统的动力响应。

图5为单个扣件失效对系统动力响应的影响。从图5（a）和图5（b）可见：由于车辆悬挂的影响，车体垂向加速度和构架垂向加速度变化并不明显，仅在扣件失效位置附近有少量增幅；另外，从图5（c）可以看出：轮轨垂向力仅在扣件失效位置附近有较大增幅，随着与扣件失效位置距离的增加，扣件失效的影响逐渐减弱；由图5（d）和图5（e）可见：扣件失效对钢轨位移和钢轨垂向加速度的影响范围较小，仅当列车到达扣件失效位置时，钢轨位移和钢轨竖向加速度有较大增幅，其峰值分别增加75.86 %和70.04%。

图5 扣件失效对车轨系统动力响应的影响

由此可见，单个扣件失效工况下车辆系统的动力响应变化很小；对于钢轨的动力响应影响较大，过大的钢轨位移和加速度会导致轨道不平顺持续恶化，增大轮轨垂向力，甚至降低钢轨的使用寿命。如图6所示。

图6 两种工况下钢轨频域响应对比图

从图6两种工况下钢轨频域响应对比图可以看出，单个扣件失效工况下钢轨垂向加速度振级在1 Hz～60 Hz频段内略大于无扣件失效时的钢轨垂向加速度振级；60 Hz以后频段两种工况的钢轨垂向加速度振级基本一致，且在125 Hz左右出现峰值。

3.2 扣件失效数量对系统振动特性的影响

扣件失效数量分别取钢轨跨中位置单个扣件失效、跨中位置扣件及其左右相邻扣件3个扣件失效，考虑扣件失效数量对系统振动特性的影响，其计算结果如图7所示。

图7 扣件失效数量对车轨系统动力响应的影响

由图7（a）可见：车体垂向加速度由无扣件失效时的0.058 m/s2增加到3 个扣件同时失效时的0.137 m/s2，增大了2.36倍，且对车体加速度的影响范围较大；由图7（b）可知：转向架垂向加速度值仅在扣件失效位置附近有较大增幅，且3 个扣件同时失效工况下的加速度幅值增大更为明显；从图7（c）可以看出：扣件失效位置处轮轨垂向力有较大增加，其他位置处增幅不明显，说明扣件失效对轮轨垂向力的影响范围较小，仅在扣件失效位置处对其有较大影响；此外，从图7（d）和图7（e）可以看出：钢轨位移随着扣件失效数量的增加而大幅增大，当3 个扣件同时失效时，其钢轨位移分别是无扣件失效和单个扣件失效时的9.67 倍和5.80 倍，与之对应的扣件失效位置处的钢轨垂向加速度峰值分别是无扣件失效、单个扣件失效时的1.65 倍和0.97 倍，达到了59.13 m/s2，虽然单个扣件失效时的钢轨加速度峰值最大，但整体钢轨加速度值与3 个扣件失效时相比，其增幅较小。因此，扣件失效数量的增加会显著增大车轨系统的动力响应，影响乘客乘车的舒适性，加剧轮轨的磨耗和钢轨扣件的失效，应及时排查检修，更换失效扣件。

从图8钢轨频域响应对比图可以看出，3种工况下钢轨垂向加速度振级在1 Hz～60 Hz 频段内相差较大；在20 Hz～50 Hz 频段内3 个扣件失效时的钢轨垂向加速度振级略小于单个扣件失效时的振级，但从整体来看，3个扣件失效时钢轨垂向加速度振级增幅最大；60 Hz以后频段内3个工况下的钢轨垂向加速度振级基本趋于一致，并在125 Hz频率左右出现峰值。

图8 3种工况下钢轨频域响应对比图

4 结语

本文在车辆-轨道耦合动力学基础上，基于ANSYS的二次开发技术，利用APDL语言提出了一种求解车辆-轨道耦合系统动力响应的方法，并以扣件失效为例，分析了扣件失效数量对车辆-轨道耦合系统动力响应的影响，得出如下结论。

（1）单个扣件失效工况下车辆系统的动力响应变化并不明显，对于轨道结构的动力响应影响较大，钢轨垂向位移和加速度值在扣件失效位置处增大明显，影响轨道结构的耐久性和安全性。

（2）车辆-轨道耦合系统的各项动力特性随着扣件失效数量的增加而发生变化，加剧轮轨的磨耗和增大临近扣件的受力，造成扣件的二次失效，影响乘客乘车的舒适性和环境振动，需及时检修。

（3）利用该方法可对不同型式的轨道结构和轨下基础进行计算分析，提出的计算方法对于研究轮轨耦合动力特性具有很好的适用性。