基于交叉迭代法的车辆-轨道非线性空间耦合振动模型及算法改进

Model and Algorithm Improvement of Vehicle-Track Nonlinear Spatial Coupling Vibration Based on Cross-Iteration Method

摘要

针对车辆-轨道非线性空间耦合系统动力学分析计算量大、计算耗时长及计算精度难以保证的问题，对车辆-轨道非线性空间耦合振动模型和交叉迭代法进行改进。运用有限元法分别建立车辆和轨道空间子系统动力学模型，在基于“迹线法”构建轮轨接触几何关系的基础上，提出“投影对点作差法”搜索轮轨空间接触点位置，并引入轮轨准弹性接触对其进行修正，精细化轮轨接触关系；结合车辆-轨道非线性空间耦合系统特点，改进基于Newmark数值积分的交叉迭代求解系统动力学方程的算法，并给出完整的数值计算步骤。通过与相关文献进行对比，验证改进模型和算法的有效性。结果表明：改进模型的分析精度更高，搜索轮轨空间接触点位置的速度更快；改进算法的计算过程更完善，步骤更清晰；改进的模型和算法在提高计算精度的同时，仍然具有较快的计算效率，且使数值编程更易实现，方便工程应用。

Abstract

In order to address the challenges posed by extensive computational requirements, prolonged computation time and reduced calculation accuracy in the dynamic analysis of vehicle-track nonlinear spatial coupling systems,the vehicle-track nonlinear spatial coupling vibration model is enhanced along with improvements made to the cross iterative method. The dynamic models for the vehicle and track spatial subsystems are separately established using the finite element method. The wheel-rail contact geometric relationship is constructed based on the "trajectory method". A "projection point difference method" is proposed to determine the position of wheel-rail spatial contact points while introducing wheel-rail quasi-elastic contact for refinement purposes in order to enhance the wheel-rail contact relationship. Considering the characteristics of vehicle-track nonlinear spatial coupling systems,the algorithm for solving system dynamic equations based on Newmark numerical integration and cross iterative methods has been improved while providing a comprehensive outline of numerical calculation steps. The effectiveness of the improved model and algorithm is verified by comparing with the calculation results of related literatures. The study shows that: this enhanced model offers superior analytical accuracy enabling faster localization of wheel-rail spatial contact points; meanwhile, the refined algorithm presents a more complete calculation process with clearer steps; both enhancements contribute to improved calculation accuracy without compromising efficiency, thus facilitating easier implementation in numerical programming applications as well as convenient engineering utilization.

Graphical abstract

关键词

车辆系统 / 轨道系统 / 非线性空间耦合 / 交叉迭代 / 轮轨接触 / 迹线法

Key words

Vehicle system / Track system / Nonlinear spatial coupling / Cross iteration / Wheel-rail contact / Trajectory method

随着高速铁路的快速发展、列车速度的不断提高以及新型车辆和新型轨道结构的大量应用，车辆与轨道间的动态相互作用更加复杂。为确保高速列车运行的安全性和平稳性，必须对车辆-轨道耦合系统在高速下的动态力学行为和服役性能进行精确分析，而车辆-轨道耦合动力学则是上述工作的基础。经过数十年的发展，车辆-轨道耦合动力学理论、模型和算法的研究日益深化^［1］，构建的系统模型精度越来越高，目前的研究方向已逐渐从垂向耦合向空间耦合过渡。相较于垂向耦合系统，空间耦合系统求解需要寻找到每一时刻轮轨接触点的位置。早期研究中，严隽耄^［2］提出了轮轨空间约束几何计算和在曲面内扫描搜索接触点的方法。王开文^［3］分析了轮轨接触几何关系特点，确定了轮轨接触点轨迹为1条在踏面上的空间曲线（迹线），使搜索轮轨接触点由二维曲面扫描降维成为一维空间曲线扫描，计算量显著降低。倪平涛等^［4-6］对迹线法进行了深入研究，给出了2种常用坐标系下轮轨接触点的正确计算式，纠正了前期研究中存在的问题，并应用赫兹非线性接触理论计算了轮轨法向接触力。高浩等^［7］提出了考虑轮对弹性的轮轨接触点计算方法。干锋^［8］提出了基于空间矢量映射求解轮轨接触点的新型算法，并将轮轨接触视为准弹性接触，对接触点进行了修正。翟婉明^［9］系统地阐述了车辆-轨道耦合动力学基本原理，给出了车辆-轨道垂向耦合振动和空间耦合振动的数学表达式和详细模型，指出采用“最小距离法”求解轮轨接触点时无须通过重复调整侧滚角和反复迭代促使左右侧轮轨最小距离相等，从而反映左右侧轮轨相互作用实时状态的不一致性。

随着研究的逐渐深入，系统模型的自由度数和复杂程度也大幅增加，进而加大了计算耗时，因此，寻求高效的求解算法一直是研究热点，且对于空间耦合复杂大系统尤为迫切。车辆-轨道耦合系统动力响应的求解方法一般有频域法和时域法2种，但频域法不适用于非线性模型动力响应的求解，若系统考虑了非线性特性，则需要在时域内求解^［10］。时域求解方法分为同步求解法和分离迭代法2种^［11］。同步求解法是先建立车辆-轨道整体耦合运动方程，使其子系统间的力协调关系自然满足，再应用逐步积分法直接求解，可避免迭代，但由于系数矩阵的时变性，每一积分步的系统矩阵都需进行更新和求逆，导致计算效率降低。分离迭代法将整个系统划分为车辆、轨道及桥梁几个非时变子系统，并通过位移协调和力平衡条件将子系统联系起来，可避免系统有效刚度矩阵每次求逆，提高求解效率，但针对计算规模较大的系统，算法存在不易收敛和迭代次数过多等问题。因此，需要结合系统的具体情况，寻找到高效的分离迭代方法，提高求解效率。雷晓燕^［12］运用有限元法建立了车辆和轨道2个子系统模型，并利用轮轨非线性接触力和位移协调条件将系统耦合，提出了基于Newmark数值计算的交叉迭代法。吴神花等^［13］引入了基于松弛法的修正因子对轮轨接触力进行修正，提高了计算收敛速度，调大了最大迭代时间步长，优化了交叉迭代法。张斌等^［14］针对车辆-轨道耦合系统运动方程联立求解过程，基于Newmark积分格式规则，对子系统方程非荷载项矩阵进行修正和求逆预处理，提出了一种将有限元法与非线性接触理论相结合的交叉迭代数值改进算法。高芒芒等^［15］采用强迫位移和强迫速度的方法处理轨下结构对钢轨系统的作用，采用施加外荷载的方式处理钢轨系统对轨下结构的作用，优化了传统的线路结构模型和传统的分离迭代算法，提出了基于强迫振动的列车-轨道-轨下结构垂向耦合动力分析方法。从以上研究中可以看出，因易于理解、便于编程及计算效率和分析精度高等优点，交叉迭代法已较为成熟地应用于求解车辆-轨道二维垂向耦合系统的动力学响应，而应用于车辆-轨道非线性空间耦合系统进行高效求解却鲜有研究。尽管雷晓燕等^［16］基于有限元法建立车辆、轨道2个空间子系统，初步探讨了应用迹线法融入交叉迭代法求解车辆-轨道非线性空间耦合系统动力学响应的算法，计算效率和分析精度有了明显提升，但是为具有更好的工程适用性，提出的模型和算法有待进一步完善。

本文运用有限元法建立车辆和轨道2个空间子系统，并基于“迹线法”原理，构建实现2个系统耦合的轮轨接触几何关系。针对搜索接触点的位置存在过程繁琐、效率不高等问题，提出“投影对点作差法”，并引入轮轨准弹性接触对接触点进行修正，精细化了轮轨接触关系。接着，改进了基于Newmark数值积分的交叉迭代求解车辆-轨道非线性空间耦合系统动力学方程的算法，给出了改进算法完整的数值计算过程。最后，通过几个数值应用说明了模型和算法的改进效果。

1 车辆-轨道非线性空间耦合系统改进模型

1.1　车辆子系统模型

车辆子系统为附有二系悬挂的整车模型，共由7个刚体部件组成，其模型如图1所示，且定义车辆子系统模型各参数变量见表1。图表中，由下标字母p区分标识各刚体部件，p取c表示车体；取t和s分别表示前后转向架；取w1，w2，w3和w4分别表示1位—4位轮对。

考虑车体和每个转向架的沉浮、横移、点头、侧滚和摇头振动，以及每个轮对的沉浮、横移、侧滚和摇头振动，车辆空间子系统共有31个自由度，详细参见文献［17］。

基于多刚体动力学原理，建立车辆子系统的动力学方程为

{\mathit{M}}_{\mathrm{u}}{\ddot{\mathit{a}}}_{\mathrm{u}}+{\mathit{C}}_{\mathrm{u}}{\dot{\mathit{a}}}_{\mathrm{u}}+{\mathit{K}}_{\mathrm{u}}{\mathit{a}}_{\mathrm{u}}={\mathit{F}}_{\mathrm{u}\mathrm{g}}+{\mathit{F}}_{\mathrm{u}}（1）

式中：{\mathit{M}}_{\mathrm{u}}，{\mathit{C}}_{\mathrm{u}}和{\mathit{K}}_{\mathrm{u}}分别为车辆子系统的质量、阻尼和刚度矩阵，其显式表达式可参见文献［18］；{\mathit{a}}_{\mathrm{u}}，{\dot{\mathit{a}}}_{\mathrm{u}}和{\ddot{\mathit{a}}}_{\mathrm{u}}分别为车辆子系统的位移、速度和加速度向量；{\mathit{F}}_{\mathrm{u}\mathrm{g}}为车辆子系统的重力向量；{\mathit{F}}_{\mathrm{u}}为车辆子系统受轮轨法向接触力作用的等效节点荷载向量，由各车轮所受的轮轨法向接触力等效变换后组成。

作用在j（j=1—4）位轮对i侧（i=1，2，其中1表示左侧，2表示右侧）车轮踏面上的轮轨法向接触力向量{\mathit{F}}_{\mathrm{u}ij}由Hertz非线性接触理论计算，为

（2）

其中，

{\partial }_{ij}=\frac{\mathrm{\Delta }{z}_{ij}}{\mathrm{c}\mathrm{o}\mathrm{s}\left[{\delta }_{ij}+{\left(-\mathrm{1}\right)}^{i+\mathrm{1}}{\theta }_{\mathrm{w}j}\right]}

\mathrm{\Delta }{z}_{ij}=\mathrm{m}\mathrm{i}\mathrm{n}\left(z_{\mathrm{r}ij}-z_{\mathrm{u}ij}-z_{\mathrm{o}ij}\right)

式中：G_{ij}为接触挠度系数，当车轮为锥形踏面时G_{ij}=\mathrm{457}{R}_{ij}^{-\mathrm{0.149}}(\mathrm{\mu }\mathrm{m}·{\mathrm{N}}^{-\mathrm{2}/\mathrm{3}})，当车轮为磨耗形踏面时G_{ij}=\mathrm{386}{R}_{ij}^{-\mathrm{0.149}}(\mathrm{\mu }\mathrm{m}·{\mathrm{N}}^{-\mathrm{2}/\mathrm{3}})，其中R_{ij}为轮轨接触点位置的车轮滚动圆半径；{\partial }_{ij}为轮轨法向接触压缩量；z_{\mathrm{r}ij}为钢轨廓形动态曲线各离散点的Z向坐标；z_{\mathrm{u}ij}为与钢轨廓形动态曲线离散点Y向坐标相等所对应车轮迹线离散点的Z向坐标；z_{\mathrm{o}ij}为车轮与钢轨无变形时（考虑轨底坡）的初始间隙；\mathrm{\Delta }{z}_{ij}为轮轨接触点位置的轮轨Z向相对位移；{\delta }_{ij}为在车轮踏面初始廓形线上对应于轮轨接触点位置的接触角；{\theta }_{\mathrm{w}j}为j位轮对侧滚角。

需要指出的是，式（2）中选用如图1（d）所示的全局坐标系中Z轴向下，采用最小距离法或最小值法确定接触点位置时应为钢轨廓形动态曲线离散点Z向坐标z_{\mathrm{r}ij}减去车轮迹线离散点Z向坐标z_{\mathrm{u}ij}。

1.2　轨道子系统模型

由于轨道的CA砂浆层和底座板具有较大的刚度，将无砟轨道子系统简化为2层轨道空间结构模型，如图2所示。基于图1（d）所示的全局坐标系，取纵向相邻扣件间的轨道结构为1个轨道空间单元，每个单元由2个空间梁单元、4个离散支承黏弹性阻尼元件、5个块体单元和连续支承黏弹性阻尼元件组成，共计28个结点，其中空间梁单元每个节点考虑X，Y和Z向的振动和绕X，Y和Z轴的转动，块体单元每个节点考虑X，Y和Z向的振动，详细参见文献［17］。

运用有限元法，建立轨道子系统的动力学方程为

{\mathit{M}}_{\mathrm{r}}{\ddot{\mathit{a}}}_{\mathrm{r}}+{\mathit{C}}_{\mathrm{r}}{\dot{\mathit{a}}}_{\mathrm{r}}+{\mathit{K}}_{\mathrm{r}}{\mathit{a}}_{\mathrm{r}}={\mathit{F}}_{\mathrm{r}\mathrm{g}}+{\mathit{F}}_{\mathrm{r}}（3）

式中：{\mathit{M}}_{\mathrm{r}}，{\mathit{C}}_{\mathrm{r}}和{\mathit{K}}_{\mathrm{r}}分别为轨道子系统的质量、阻尼和刚度矩阵；{\mathit{a}}_{\mathrm{r}}，{\dot{\mathit{a}}}_{\mathrm{r}}和{\ddot{\mathit{a}}}_{\mathrm{r}}分别为轨道子系统的位移、速度和加速度向量；{\mathit{F}}_{\mathrm{r}\mathrm{g}}为轨道子系统重力向量；{\mathit{F}}_{\mathrm{r}}为轨道子系统受轮轨法向接触力作用的等效节点荷载向量，由钢轨在j位轮对i侧车轮接触点位置承受的轮轨法向接触力{\mathit{F}}_{\mathrm{r}ij}等效变换后组成，{\mathit{F}}_{\mathrm{r}ij}与{\mathit{F}}_{\mathrm{u}ij}大小相等、方向相反。

1.3　轮轨接触关系

轮轨接触关系是车辆子系统和轨道子系统的连接纽带^［19］，是耦合系统求解的关键，进行轮轨接触几何计算通常采用迹线法^［4］。基于迹线法原理，约定左侧和右侧滚动圆横坐标d_{\mathrm{u}ij}考虑正负，接触角{\delta }_{ij}均为正值，推导生成某时刻图1（d）所示全局坐标系下车轮空间接触点的坐标（X_{\mathrm{u}ij}，Y_{\mathrm{u}ij}，Z_{\mathrm{u}ij}）为

\left\{\begin{array}{c}{X}_{\mathrm{u}ij}=X_{\mathrm{o}ij}+{(-\mathrm{1})}^{i+\mathrm{1}}{l}_{xj}{R}_{ij}\mathrm{t}\mathrm{a}\mathrm{n}{\delta }_{ij}+X_{\mathrm{w}j}\\ \begin{array}{l}{Y}_{\mathrm{u}ij}=Y_{\mathrm{o}ij}+{(-\mathrm{1})}^i\frac{R_{ij}}{\mathrm{1}-{l_{xj}}^{\mathrm{2}}}[{l_{xj}}^{\mathrm{2}}{l}_{yj}\mathrm{t}\mathrm{a}\mathrm{n}{\delta }_{ij}+\\ {(-\mathrm{1})}^{i+\mathrm{1}}{l}_{zj}{m}_{ij}]+Y_{\mathrm{w}j}\end{array}\\ \begin{array}{l}{Z}_{\mathrm{u}ij}=Z_{\mathrm{o}ij}+{(-\mathrm{1})}^i\frac{R_{ij}}{\mathrm{1}-{l_{xj}}^{\mathrm{2}}}[{l_{xj}}^{\mathrm{2}}{l}_{zj}\mathrm{t}\mathrm{a}\mathrm{n}{\delta }_{ij}+\\ {(-\mathrm{1})}^i{l}_{yj}{m}_{ij}]+Z_{\mathrm{w}j}\end{array}\end{array}\right.（4）

其中，

X_{\mathrm{o}ij}=d_{\mathrm{u}ij}{l}_{xj}

Y_{\mathrm{o}ij}=d_{\mathrm{u}ij}{l}_{yj}

Z_{\mathrm{o}ij}=d_{\mathrm{u}ij}{l}_{zj}

l_{xj}=-\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{w}j}\mathrm{s}\mathrm{i}\mathrm{n}{\psi }_{\mathrm{w}j}

l_{yj}=\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{w}j}\mathrm{c}\mathrm{o}\mathrm{s}{\psi }_{\mathrm{w}j}

l_{zj}=\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{w}j}

m_{ij}=\sqrt[]{\mathrm{1}-{l_{xj}}^{\mathrm{2}}\left(\mathrm{1}+\mathrm{t}\mathrm{a}{\mathrm{n}}^{\mathrm{2}}{\delta }_{ij}\right)}

式中：X_{\mathrm{o}ij}，Y_{\mathrm{o}ij}，Z_{\mathrm{o}ij}为j位轮对i侧车轮的滚动圆圆心在轮对平移坐标系中的坐标；l_{xj}，l_{yj}，l_{zj}为j位轮对轴线在轮对平移坐标系中的方向余弦；{\theta }_{\mathrm{w}j}，{\psi }_{\mathrm{w}j}和X_{\mathrm{w}j}，Y_{\mathrm{w}j}，Z_{\mathrm{w}j}分别为j位轮对质心在全局坐标系中的侧滚角、摇头角和X，Y和Z向位移；d_{\mathrm{u}ij}为j位轮对i侧车轮踏面滚动圆圆心在轮对本体坐标系中的横向坐标。

将某时刻轮对位移参数和车轮廓形离散点坐标数据代入式（4），即可生成此时刻轮轨接触点在全局坐标系中的“轨迹线”。结合轮轨接触区域的实际情况，计算时只须截用车轮踏面廓形的黄实线部分（如图3所示），这样既不影响后续计算结果，又避免了接触角{\delta }_{ij}正负值混淆问题，还适当降低了计算量。

图1（d）所示全局坐标系的原点不在钢轨截面中心点的延长线上，因此，需先获取钢轨截面中心与图1（d）全局坐标系原点重合时的轨头廓形离散点坐标数据数组。依据钢轨位移响应、轨道不平顺和轨道结构等参数（轨距为考虑轨底坡时钢轨轨顶面下16 mm范围内两股工作边之间的最小距离，而轨底坡应以钢轨截面中心为旋转点），先旋转后平移得到在全局坐标系中钢轨轨头某时刻的廓形动态曲线，其表达式为

\left\{\begin{array}{c}\begin{array}{l}{Y}_{\mathrm{r}ij}=Y_{\mathrm{r}\mathrm{q}i}\mathrm{c}\mathrm{o}\mathrm{s}[{\theta }_{\mathrm{r}\mathrm{w}ij}+{(-\mathrm{1})}^{i+\mathrm{1}}{\theta }_{\mathrm{r}\mathrm{o}}]-\\ Z_{\mathrm{r}\mathrm{q}i}\mathrm{s}\mathrm{i}\mathrm{n}[{\theta }_{\mathrm{r}\mathrm{w}ij}+{(-\mathrm{1})}^{i+\mathrm{1}}{\theta }_{\mathrm{r}\mathrm{o}}]+\\ Y_{\mathrm{r}\mathrm{o}i}+Y_{\mathrm{r}\mathrm{w}ij}+{\eta }_{\mathrm{r}ij}\end{array}\\ \begin{array}{l}{Z}_{\mathrm{r}ij}=Y_{\mathrm{r}\mathrm{q}i}\mathrm{s}\mathrm{i}\mathrm{n}[{\theta }_{\mathrm{r}\mathrm{w}ij}+{(-\mathrm{1})}^{i+\mathrm{1}}{\theta }_{\mathrm{r}\mathrm{o}}]+\\ Z_{\mathrm{r}\mathrm{q}i}\mathrm{c}\mathrm{o}\mathrm{s}[{\theta }_{\mathrm{r}\mathrm{w}ij}+{(-\mathrm{1})}^{i+\mathrm{1}}{\theta }_{\mathrm{r}\mathrm{o}}]+\\ Z_{\mathrm{r}\mathrm{o}i}+Z_{\mathrm{r}\mathrm{w}ij}+{\zeta }_{\mathrm{r}ij}\end{array}\end{array}\right.（5）

式中：Y_{\mathrm{r}ij}和Z_{\mathrm{r}ij}分别为j位轮对i侧轮轨接触点位置截面的钢轨轨头廓形动态曲线离散点Y向和Z向坐标数组；Y_{\mathrm{r}\mathrm{q}i}和Z_{\mathrm{r}\mathrm{q}i}分别为在钢轨局部坐标系下i侧钢轨无变形时轨头廓形离散点的Y和Z向坐标数组；{\theta }_{\mathrm{r}\mathrm{w}ij}和Y_{\mathrm{r}\mathrm{w}ij}，Z_{\mathrm{r}\mathrm{w}ij}分别为j位轮对i侧轮轨接触点位置的钢轨截面中心的侧滚角及Y向和Z向位移；{\mathit{\theta }}_{\mathrm{r}\mathrm{o}}为钢轨轨底坡；{\eta }_{\mathrm{r}ij}和{\zeta }_{\mathrm{r}ij}分别为j位轮对i侧轮轨接触点位置钢轨的横向和高低不平顺；Y_{\mathrm{r}\mathrm{o}i}和Z_{\mathrm{r}\mathrm{o}i}分别为i侧钢轨截面中心在全局坐标系中的Y向和Z向坐标。

在构建了上述轮轨空间接触几何关系后，以往研究^［4，16］以迹线离散点数量为次数，通过反复循环扫描比较各离散点对应的轮轨垂向距离，最终通过“垂向最小距离”确定轮轨接触点位置。该办法搜索接触点的过程繁琐，而且迹线离散密度越大，循环搜索次数越多，计算越耗时。

1.4　基于投影对点作差法的轮轨接触点计算及修正

1.4.1　算法提出

针对上述问题，提出一种高效搜索轮轨接触点位置的新办法——投影对点作差法。首先，截取车轮空间接触点迹线在钢轨廓形动态曲线Y坐标范围内的线段，并以此线段离散点的Y坐标重新插值离散钢轨廓形动态曲线，生成与迹线截取线段Y坐标相等、Z坐标一一对应的钢轨廓形动态曲线新离散点坐标数组；其次，计算2个数组的Z向坐标差并寻找最小差值，该最小差值即为轮轨接触点处的轮轨Z向相对位移，而该最小差值对应的离散点即为轮轨接触点；然后，利用接触点的Y坐标，找到接触点在原完整车轮迹线离散点数组内的位置。

因为图形的平移和旋转变换不改变离散点坐标在数组内的位置，所以可按接触点在数组内的位置直接调用所需要的车轮参数。另外，钢轨参数可以通过接触点的Y坐标插值获得，也可以近似调用钢轨廓形动态曲线初始离散点中与实际接触点位置最为接近点的参数。需要指出的是，由于钢轨廓形曲线比较平滑，当钢轨初始离散密度较大时，采用上述接近点参数产生的误差很小。“投影对点作差法”通过构建车轮和钢轨在接触区域内的相同大小坐标数组后再计算2个数组的Z向坐标差，避免了沿空间迹线逐点扫描，可大幅提高搜索接触点位置的效率。

1.4.2　算法验证

以LM_A型车轮和60 kg · m^-1级钢轨为例，采用“投影对点作差法”进行搜索接触点位置计算。车轮踏面和钢轨轨头廓形由分段曲线函数拟合生成，具体函数表达式可参见文献［20-22］。曲线函数可对横坐标按某一密度离散，生成曲线离散点及其坐标数据数组，以此生成的钢轨轨头廓形（右侧）和车轮踏面型面（右侧）分别如图4和图5所示。接着，在图1（d）所示的全局坐标系下依据车辆-轨道非线性空间耦合系统结构参数，生成左右车轮踏面和钢轨轨头廓形离散点坐标数据数组。

通过“投影对点作差法”计算得到的轮轨刚性接触点位置分布和接触角与文献［23］的对比结果分别如图6和图7所示。由图6和图7可知：二者的计算结果吻合良好，接触点位置分布和接触角变化趋势基本一致，验证了本文算法的有效性；局部差异可能是车轮和钢轨廓型曲线存在拟合误差和曲线离散密度不同引起的。

从图6（a）还可知：将轮轨接触考虑成刚性接触，则接触点位置不变形，轮轨顶部接触区域会有较大的位置跳跃。相比较动力学主流商业软件SIMPACK，其轮轨关系模块考虑了轮轨接触点处的弹性变形，但未考虑轮对横移时的侧滚振动，无法反映轮轨真实的接触状态^［24］。为不失计算的准确性，本文在提出能同时兼顾轮对横移和侧滚运动的“投影对点作差法”基础上，视轮轨接触为准弹性接触，对计算得到的刚性接触点位置进行修正。准弹性修正法可参见文献［8］，对数应变系数\epsilon取4.5×10^-5。

准弹性修正后轮轨接触点位置分布对比结果如图8所示。从图8可知：本文计算结果与文献［8］吻合良好，证实了本文办法的可行性；此外，与图6（a）相比，图8（a）显示的轮轨接触关系平滑、均匀、连续，原轮轨顶部接触区域存在的位置跳跃现象得到了明显改善；值得一提的是，准弹性修正法能有效地提高“投影对点作差法”的计算精度，但对其计算效率几乎没有负影响，二者匹配甚佳。

2 车辆-轨道耦合系统动力学求解的交叉迭代改进算法

2.1　计算步骤

由于车辆-轨道非线性空间耦合系统动力学方程是时变的，采用Newmark数值积分法分别求解车辆子系统和轨道子系统动力学方程。在每一计算时间步中，对车辆子系统和轨道子系统动力学方程交叉迭代求解。因此，需将车辆子系统和轨道子系统动力学方程按Newmark积分格式预处理，得到由t-\mathrm{\Delta }t时刻的动力学响应表示t时刻第e次迭代动力学响应的表达式如下。

车辆子系统位移表达式为

\begin{array}{l}({\mathit{K}}_{\mathrm{u}}+c_{\mathrm{0}}{\mathit{M}}_{\mathrm{u}}+c_{\mathrm{1}}{\mathit{C}}_{\mathrm{u}}){\mathit{a}}_{\mathrm{u}e}(t)=\\ {\mathit{F}}_{\mathrm{u}\mathrm{g}}+{\mathit{F}}_{\mathrm{u}e}(t)+{\mathit{M}}_{\mathrm{u}}(c_{\mathrm{0}}{\mathit{a}}_{\mathrm{u}}(t-\mathrm{\Delta }t)+\\ c_{\mathrm{2}}{\dot{\mathit{a}}}_{\mathrm{u}}(t-\mathrm{\Delta }t)+c_{\mathrm{3}}{\ddot{\mathit{a}}}_{\mathrm{u}}(t-\mathrm{\Delta }t))+\\ {\mathit{C}}_{\mathrm{u}}(c_{\mathrm{1}}{\mathit{a}}_{\mathrm{u}}(t-\mathrm{\Delta }t)+c_{\mathrm{4}}{\dot{\mathit{a}}}_{\mathrm{u}}(t-\mathrm{\Delta }t)+\\ c_{\mathrm{5}}{\ddot{\mathit{a}}}_{\mathrm{u}}(t-\mathrm{\Delta }t))\end{array}（6）

车辆子系统加速度表达式为

\begin{array}{l}{\ddot{\mathit{a}}}_{\mathrm{u}e}(t)=c_{\mathrm{0}}({\mathit{a}}_{\mathrm{u}e}(t)-{\mathit{a}}_{\mathrm{u}}(t-\mathrm{\Delta }t))-\\ c_{\mathrm{2}}{\dot{\mathit{a}}}_{\mathrm{u}}(t-\mathrm{\Delta }t)-c_{\mathrm{3}}{\ddot{\mathit{a}}}_{\mathrm{u}}(t-\mathrm{\Delta }t)\end{array}（7）

车辆子系统速度表达式为

{\dot{\mathit{a}}}_{\mathrm{u}e}(t)={\dot{\mathit{a}}}_{\mathrm{u}}(t-\mathrm{\Delta }t)+c_{\mathrm{6}}{\ddot{\mathit{a}}}_{\mathrm{u}}(t-\mathrm{\Delta }t)+c_{\mathrm{7}}{\ddot{\mathit{a}}}_{\mathrm{u}e}(t)（8）

轨道子系统位移表达式为

\begin{array}{l}({\mathit{K}}_{\mathrm{r}}+c_{\mathrm{0}}{\mathit{M}}_{\mathrm{r}}+c_{\mathrm{1}}{\mathit{C}}_{\mathrm{r}}){\mathit{a}}_{\mathrm{r}e}(t)=\\ {\mathit{F}}_{\mathrm{r}\mathrm{g}}+{\mathit{F}}_{\mathrm{r}e}(t)+{\mathit{M}}_{\mathrm{r}}(c_{\mathrm{0}}{\mathit{a}}_{\mathrm{r}}(t-\mathrm{\Delta }t)+\\ c_{\mathrm{2}}{\dot{\mathit{a}}}_{\mathrm{r}}(t-\mathrm{\Delta }t)+c_{\mathrm{3}}{\ddot{\mathit{a}}}_{\mathrm{r}}(t-\mathrm{\Delta }t))+\\ {\mathit{C}}_{\mathrm{r}}(c_{\mathrm{1}}{\mathit{a}}_{\mathrm{r}}(t-\mathrm{\Delta }t)+c_{\mathrm{4}}{\dot{\mathit{a}}}_{\mathrm{r}}(t-\mathrm{\Delta }t)+\\ c_{\mathrm{5}}{\ddot{\mathit{a}}}_{\mathrm{r}}(t-\mathrm{\Delta }t))\end{array}（9）

轨道子系统加速度表达式为

\begin{array}{l}{\ddot{\mathit{a}}}_{\mathrm{r}e}(t)=c_{\mathrm{0}}({\mathit{a}}_{\mathrm{r}e}(t)-{\mathit{a}}_{\mathrm{r}}(t-\mathrm{\Delta }t))-\\ c_{\mathrm{2}}{\dot{\mathit{a}}}_{\mathrm{r}}(t-\mathrm{\Delta }t)-c_{\mathrm{3}}{\ddot{\mathit{a}}}_{\mathrm{r}}(t-\mathrm{\Delta }t)\end{array}（10）

轨道子系统速度表达式为

{\dot{\mathit{a}}}_{\mathrm{r}e}(t)={\dot{\mathit{a}}}_{\mathrm{r}}(t-\mathrm{\Delta }t)+c_{\mathrm{6}}{\ddot{\mathit{a}}}_{\mathrm{r}}(t-\mathrm{\Delta }t)+c_{\mathrm{7}}{\ddot{\mathit{a}}}_{\mathrm{r}e}(t)（11）

式中：c_{\mathrm{m}}\left(m=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4},\mathrm{5},\mathrm{6},\mathrm{7}\right)为Newmark数值计算法的积分常数^［25］；{\mathit{a}}_{\mathrm{u}e}(t)，{\dot{\mathit{a}}}_{\mathrm{u}e}(t)，{\ddot{\mathit{a}}}_{\mathrm{u}e}(t)和{\mathit{F}}_{\mathrm{u}e}(t)分别为t时刻第e次迭代时车辆子系统的位移、速度、加速度向量和受轮轨法向接触力作用的等效节点荷载向量；{\mathit{a}}_{\mathrm{r}e}(t)，{\dot{\mathit{a}}}_{\mathrm{r}e}(t)，{\ddot{\mathit{a}}}_{\mathrm{r}e}(t)和{\mathit{F}}_{\mathrm{r}e}(t)分别为t时刻第e次迭代时轨道子系统的位移、速度、加速度向量和受轮轨法向接触力作用的等效节点荷载向量。

由于车辆子系统和轨道子系统方程式中的有效刚度矩阵为非时变矩阵，可先求逆预处理，后续计算直接调用即可，避免了有效刚度矩阵多次求逆，提高了计算效率。在文献［12］给出的算法步骤中，e循环的每一步，需进行2次轮轨力求解，即第1次求解得到轮轨力后计算轨道子系统振动响应，然后以此轨道新响应值，再一次求解轮轨力，并将此轮轨力作用于车辆子系统。该过程在求解垂向耦合振动方程时，可降低e循环的最大迭代次数，提高计算效率。对于空间耦合系统，求解轮轨力须先寻找到轮轨接触点的位置，且求解2次轮轨力须寻找2次轮轨接触点位置，这样e循环的每一步中，存在2个接触点，只有当这2个接触点位置接近重合时，计算才会收敛，这势必导致迭代次数增加，影响计算效率。因此，本文结合轮轨空间耦合关系的具体情况，对原有算法作了进一步改进。现给出应用交叉迭代改进算法的主要计算步骤如下。

2.1.1　初始计算

（1）初始计算（第零时间步）时，假设车辆和轨道子系统结构初始状态为静止，即位移、速度和加速度为零，则轮轨接触点位置的钢轨截面中心的位移也自然为零。

（2）由于Newmark数值积分法需已知t-\mathrm{\Delta }t时刻的动力学响应才能求解t时刻动力学响应，因此，需要对第一时间步中e循环的首次迭代进行特殊处理。首次迭代时，轨道子系统为初始静止状态，车辆子系统则考虑受自重和轮轨接触力作用，且轮轨接触刚度作线性化处理，则车辆子系统动力响应表达式为

{\mathit{M}}_{\mathrm{u}}{\ddot{\mathit{a}}}_{\mathrm{u}}+{\mathit{C}}_{\mathrm{u}}{\dot{\mathit{a}}}_{\mathrm{u}}+({\mathit{K}}_{\mathrm{u}}+{\mathit{K}}_{\mathrm{w}}){\mathit{a}}_{\mathrm{u}}={\mathit{F}}_{\mathrm{u}\mathrm{g}}（12）

其中，

k_{\mathrm{w}}=\frac{\mathrm{3}}{\mathrm{2}{\mathrm{G}}_{\mathrm{0}}}{P_{\mathrm{0}}}^{\frac{\mathrm{1}}{\mathrm{3}}}

式中：{\mathit{K}}_{\mathrm{w}}为车辆子系统轮轨线性化接触刚度矩阵，由元素k_{\mathrm{w}}组成；G_{\mathrm{0}}为车轮名义滚动半径位置的接触挠度系数；P_{\mathrm{0}}为单轮静荷载。

（3）根据步骤（1）和步骤（2），对车辆子系统进行第一时间步中的首次迭代计算，其Newmark数值积分方程为

\begin{array}{l}({\mathit{K}}_{\mathrm{u}}+{\mathit{K}}_{\mathrm{w}}+c_{\mathrm{0}}{\mathit{M}}_{\mathrm{u}}+c_{\mathrm{1}}{\mathit{C}}_{\mathrm{u}}){\mathit{a}}_{\mathrm{u}\mathrm{1}}(\mathrm{\Delta }t)=\\ {\mathit{F}}_{\mathrm{u}\mathrm{g}}+{\mathit{F}}_{\mathrm{u}\mathrm{1}}(\mathrm{\Delta }t)+{\mathit{M}}_{\mathrm{u}}(c_{\mathrm{0}}{\mathit{a}}_{\mathrm{u}}(\mathrm{0})+\\ c_{\mathrm{2}}{\dot{\mathit{a}}}_{\mathrm{u}}(\mathrm{0})+c_{\mathrm{3}}{\ddot{\mathit{a}}}_{\mathrm{u}}(\mathrm{0}))+{\mathit{C}}_{\mathrm{u}}(c_{\mathrm{1}}{\mathit{a}}_{\mathrm{u}}(\mathrm{0})+\\ c_{\mathrm{4}}{\dot{\mathit{a}}}_{\mathrm{u}}(\mathrm{0})+c_{\mathrm{5}}{\ddot{\mathit{a}}}_{\mathrm{u}}(\mathrm{0}))\end{array}（13）

式中：{\mathit{a}}_{\mathrm{u}\mathrm{1}}(\mathrm{\Delta }t)，{\mathit{F}}_{\mathrm{u}\mathrm{1}}(\mathrm{\Delta }t)分别为第一时间步中首次迭代时车辆子系统的位移和受轮轨法向接触力作用的等效节点荷载向量。

将车辆子系统系统响应初始值代入式（13）后得

{\mathit{a}}_{\mathrm{u}\mathrm{1}}(\mathrm{\Delta }t)=({\mathit{K}}_{\mathrm{u}}+{\mathit{K}}_{\mathrm{w}}+c_{\mathrm{0}}{\mathit{M}}_{\mathrm{u}}+c_{\mathrm{1}}{\mathit{C}}_{\mathrm{u}}{)}^{-\mathrm{1}}{\mathit{F}}_{\mathrm{u}\mathrm{g}}（14）

启步后便可以进入时间步循环中，并继续第一时间步中剩余e循环迭代直至收敛，所得到的迭代稳定值即为第一时间步的最终响应量，该确定性响应量将作为下一时间步内e循环迭代的初始几何和力学边界条件。如此以时间步中内套e循环迭代的形式重复，直至计算结束。

2.1.2　时间步长循环

设已计算至第n时间步，且已进行了e-\mathrm{1}次迭代，现考察第e次迭代的计算过程，具体步骤如下。

（1）结合轮轨接触关系理论，调用在第n时间步第e-\mathrm{1}次迭代时轮轨接触点处的车轮和钢轨截面中心的位移{\mathit{a}}_{\mathrm{u}ije-\mathrm{1}}(n∆t)和{\mathit{a}}_{\mathrm{r}ije-\mathrm{1}}(n∆t)以及不平顺等参数，便可确定轮轨接触点位置，计算得到在第n时间步第e次迭代时的轮轨法向压缩量{\partial }_{ije}和法向接触力{\mathit{F}}_{\mathrm{u}ije}(n∆t)。

（2）运用松弛法对轮轨法向接触力进行修正，令

{\mathit{F}}_{\mathrm{u}ije}(n∆t)={\mathit{F}}_{\mathrm{u}ije-\mathrm{1}}(n∆t)+\mu ({\mathit{F}}_{\mathrm{u}ije}(n∆t)-{\mathit{F}}_{\mathrm{u}ije-\mathrm{1}}(n∆t))（15）

{\mathit{F}}_{\mathrm{r}ije}(n∆t)=-{\mathit{F}}_{\mathrm{u}ije}(n∆t)（16）

式中：\mu为修正因子，视时间步长取值（通常取0.3~0.5）；{\mathit{F}}_{\mathrm{u}ije}(n∆t)，{\mathit{F}}_{\mathrm{r}ije}(n∆t)分别为在第n时间步第e次迭代时的车辆、轨道子系统承受的轮轨法向接触力；{\mathit{F}}_{\mathrm{u}ije-\mathrm{1}}(n∆t)为在第n时间步第e-1次迭代时的车辆子系统承受的轮轨法向接触力。

（3）将各轮轨法向接触力{\mathit{F}}_{\mathrm{u}ije}(n∆t)，{\mathit{F}}_{\mathrm{r}ije}(n∆t)变换到全局坐标系中，并根据轮对和钢轨的运动方程，分别等效生成在第n时间步第e次迭代时施加于车辆子系统和轨道子系统的等效节点荷载向量{\mathit{F}}_{\mathrm{u}e}(n∆t)和{\mathit{F}}_{\mathrm{r}e}(n∆t)。再运用Newmark数值积分法，计算得到车辆子系统在第n时间步第e次迭代时的位移、速度和加速度响应{\mathit{a}}_{\mathrm{u}e}(n∆t)，{\dot{\mathit{a}}}_{\mathrm{u}e}(n∆t)，{\ddot{\mathit{a}}}_{\mathrm{u}e}(n∆t)，以及轨道子系统的位移、速度和加速度响应{\mathit{a}}_{\mathrm{r}e}(n∆t)，{\dot{\mathit{a}}}_{\mathrm{r}e}(n∆t)，{\ddot{\mathit{a}}}_{\mathrm{r}e}(n∆t)。

（4）计算轨道位移差值，为

∆{\mathit{a}}_{\mathrm{r}e}(n∆t)={\mathit{a}}_{\mathrm{r}e}(n∆t)-{\mathit{a}}_{\mathrm{r}e-\mathrm{1}}(n∆t)（17）

式中：{\mathit{a}}_{\mathrm{r}e-\mathrm{1}}(n∆t)为在第n时间步第e-\mathrm{1}次迭代时的轨道子系统位移。

（5）对轨道位移进行收敛性判别，收敛准则可参见文献［13］。

如果不收敛，直接转至步骤（1），再令e=e+\mathrm{1}，进入下一迭代步循环。如果收敛，令

{\mathit{a}}_{\mathrm{r}\mathrm{1}}((n+\mathrm{1})∆t)={\mathit{a}}_{\mathrm{r}e}(n∆t)（18）

{\dot{\mathit{a}}}_{\mathrm{r}\mathrm{1}}((n+\mathrm{1})∆t)={\dot{\mathit{a}}}_{\mathrm{r}e}(n∆t)（19）

{\ddot{\mathit{a}}}_{\mathrm{r}\mathrm{1}}((n+\mathrm{1})∆t)={\ddot{\mathit{a}}}_{\mathrm{r}e}(n∆t)（20）

{\mathit{a}}_{\mathrm{u}\mathrm{1}}((n+\mathrm{1})∆t)={\mathit{a}}_{\mathrm{u}e}(n∆t)（21）

{\dot{\mathit{a}}}_{\mathrm{u}\mathrm{1}}((n+\mathrm{1})∆t)={\dot{\mathit{a}}}_{\mathrm{u}e}(n∆t)（22）

{\ddot{\mathit{a}}}_{\mathrm{u}\mathrm{1}}((n+\mathrm{1})∆t)={\ddot{\mathit{a}}}_{\mathrm{u}e}(n∆t)（23）

{\mathit{a}}_{\mathrm{r}ij\mathrm{1}}((n+\mathrm{1})∆t)={\mathit{a}}_{\mathrm{r}ije}(n∆t)（24）

式中：{\mathit{a}}_{\mathrm{u}\mathrm{1}}((n+\mathrm{1})∆t)，{\dot{\mathit{a}}}_{\mathrm{u}\mathrm{1}}((n+\mathrm{1})∆t)和{\ddot{\mathit{a}}}_{\mathrm{u}\mathrm{1}}((n+\mathrm{1})∆t)分别为车辆子系统在第n+1时间步第1次迭代时的位移、速度和加速度；{\mathit{a}}_{\mathrm{r}\mathrm{1}}((n+\mathrm{1})∆t)，{\dot{\mathit{a}}}_{\mathrm{r}\mathrm{1}}((n+\mathrm{1})∆t)和{\ddot{\mathit{a}}}_{\mathrm{r}\mathrm{1}}((n+\mathrm{1})∆t)分别为轨道子系统在第n+1时间步第1次迭代时的位移、速度和加速度；{\mathit{a}}_{\mathrm{r}ij\mathrm{1}}((n+\mathrm{1})∆t)为在第n+1时间步第1次迭代时轮轨接触点处的钢轨截面中心的位移；{\mathit{a}}_{\mathrm{r}ije}(n∆t)为在第n时间步第e次迭代时轮轨接触点处的钢轨截面中心的位移；

即将第n时间步的2个子系统节点最终响应量作为第n+1时间步内e循环迭代的初始几何和力学边界条件，然后进入第n+1时间步计算。对比文献［12］的算法计算流程，本文给出的上述算法步骤中，在e循环的每一步只进行了1次轮轨力计算，然后将轮轨力以力相互作用的形式分别等效加载到车辆和轨道子系统中进行系统响应计算。这样可降低e循环的迭代次数，加速收敛。同时，可减少搜索接触点的次数，大幅降低计算耗时，提高求解效率。此外，根据有限元法原理可知，轮轨接触点位置的钢轨截面中心位移由钢轨节点位移通过插值函数求得。因此，本文在e循环迭代中，增设了式（24），这样避免了轮轨接触点位置的钢轨截面中心位移在每一时间步的迭代初值取零，能显著地提升计算效率，且使算法更加符合理论原理。

2.2　算例验证

依据德国高速轨道不平顺谱，模拟生成线路空间不平顺随机样本，并以此作为激励，计算单节车辆以200 km · h^-1速度运行工况时钢轨截面中心的振动响应，并与文献［26］的结果对比分别如图9和图10所示。图中：红色圆圈处为最大值，红色三角形处为最小值。由图9和图10可知：本文计算的位移和加速度响应均在通常范围内，响应波形也符合物理概念，说明改进模型可以较好地反映板式轨道振动特性；此外，与文献［26］相比，二者在振动趋势和波形上具有较好的一致性，计算响应幅值也与文献［26］列车的单节车辆振动平均幅值接近，验证了本文动力学耦合模型和算法的有效性；局部存在的差异可能是由模型简化和参数取值差异、轨道不平顺和动车组编组数不同引起的。

3 改进效果对比

为说明模型和算法的改进效果，与文献［16］的计算精度和计算效率进行对比。另外，轨道不平顺的随机性可能导致计算结果的不确定性，无法准确辨识二者差别，因此，下面对比分析无特殊说明时均假设轨道结构为理想平顺工况。

计算时线路长度为187.5 m，列车速度为300 km · h^-1。

3.1　计算精度对比

本文与文献［16］相同算例的钢轨截面中心垂向位移时程曲线对比计算结果如图11所示。从图11可知：平顺工况下二者计算得到的钢轨垂向位移时程曲线高度吻合，但幅值存在差异，这主要是本文对轮轨接触点位置进行了准弹性修正引起的，表明本文计算结果更准确；另外，二者计算的幅值差异并不大，这是由于未磨耗的车轮踏面和钢轨廓形曲线比较平滑，进行准弹性修正前、后计算的2个接触点位置的参数值接近所致。

3.2　计算效率对比

3.2.1　不同离散密度工况下

钢轨和车轮的廓形离散密度直接关系轮轨接触点的位置，影响轮轨力计算精度，离散密度越大，计算得到的轮轨力精度越高。

统计同台电脑多次计算同一离散密度工况的平均耗时，得到离散密度对计算耗时的影响见表2。表中：算法A为文献［16］采用的迹线法+交叉迭代法；算法B为投影对点作差法+交叉迭代法；算法C为投影对点作差法+改进的交叉迭代法。由表2可知：当钢轨和车轮的廓形离散从0.500 mm增密至0.200 mm时，文献［16］的计算耗时增加了39.458 h，而算法B和算法C的计算耗时仅增加了0.040和0.069 h；当廓形离散增密至0.001 mm时，算法B和算法C的计算耗时也仅增加了0.747和0.307 h，表明“投影对点作差法”对离散密度不敏感。

不同算法计算耗时对比如图12所示。从图12可知：算法B和算法C的计算耗时远低于算法A，且算法C的耗时又低于算法B；随着轮轨廓形离散密度的增大，算法A的耗时增加，增幅速率变大，尽管算法B和算法C的耗时也有所增加，但是增幅很小，由此可知，“投影对点作差法”的计算效率很高；从以上意义上讲，不论“迹线扫描”还是早期的“曲面扫描”，采用“投影对点作差法”均可大幅提升搜索接触点位置的效率，并且预估二者的计算耗时应该比较接近；此外，对比算法B和算法C可知，改进的交叉迭代法更能发挥算法优势，显著地提升计算效率。因此，本文提出的算法C能够在提高计算精度的同时仍然具有较快的计算效率，实现了计算精度和计算效率的最佳平衡。

3.2.2　不同收敛精度工况下

以廓形离散密度为0.500 mm时考察收敛精度对计算耗时和最大迭代次数的影响，结果见表3。从表3可知：收敛精度在10^-5~10^-8 m常规取值范围内，随着收敛精度取值的减小，算法A，算法B和算法C的最大迭代次数增多，计算耗时也会增加，但增幅不大，这是因为收敛精度取值主要在计算未稳定的初始阶段起作用，主要影响初始阶段的计算耗时，但在中、后期计算响应达到稳定后，收敛精度取值对计算耗时的影响很小，且中、后期计算为整个计算的主要部分，占比较大，致使收敛精度在不同取值工况下的计算耗时变化不大。

4 结论

（1）考虑车轮踏面廓形和钢轨轨头廓形，联合应用迹线法和准弹性修正法，构建更为精细的轮轨空间接触几何关系，并提出通过构建车轮和钢轨在接触区域内的相同大小坐标数组后再计算2个数组的Z向坐标差的搜索轮轨接触点位置的新方法——“投影对点作差法”，大幅提升了搜索轮轨接触点的效率和精确度。

（2）改进了基于Newmark数值积分的交叉迭代求解车辆-轨道非线性空间耦合系统动力学方程的算法，实现了同步进行轮轨接触点搜索与车辆-轨道非线性空间耦合系统方程求解，提高了交叉迭代收敛速度和数值分析效率。

（3）增加钢轨、车轮廓形离散密度和交叉迭代收敛精度能够提高计算精度。在文章算例中，当廓形离散密度取值从0.500 mm变成0.200 mm时，相同计算资源条件下，算法C的计算耗时从0.438 h增加到0.507 h，耗时增加极小。而文献［16］的计算耗时从16.575 h增加到了56.033 h，计算效率显著低于算法C。因此，本文提出的改进模型和算法能够在提高计算精度的同时仍然具有较快的计算效率，实现了计算精度和计算效率的最佳平衡。

（4）本文结合空间耦合系统的具体情况，对基于交叉迭代法的车辆-轨道非线性空间耦合振动模型及算法作了进一步改进，并给出了改进算法更完善的数值计算过程和更清晰的计算步骤，使数值编程更容易实现，方便工程应用。

参考文献

[1] 翟婉明.车辆-轨道耦合动力学研究的新进展［J］.中国铁道科学，2002，23（2）：1-14.

[2] ZHAI Wanming. New Advance in Vehicle-Track Coupling Dynamics ［J］. China Railway Science， 2002， 23 （2）： 1-14. in Chinese

[3] 严隽耄.具有任意轮廓形状的轮轨空间几何约束的研究［J］.西南交通大学学报，1983，30（3）：40-47.

[4] YAN Juanmao. A Study on Wheel/Rail Spatial Geometric Constraints for Arbitrary Wheel and Rail Profiles ［J］. Journal of Southwest Jiaotong University， 1983， 30 （3）： 40-47. in Chinese

[5] 王开文.车轮接触点迹线及轮轨接触几何参数的计算［J］.西南交通大学学报，1984，31（1）：89-98.

[6] WANG Kaiwen. The Track of Wheel Contact Points and the Calculation of Wheel/Rail Geometric Contact Parameters ［J］. Journal of Southwest Jiaotong University， 1984， 31 （1）： 89-98. in Chinese

[7] 倪平涛，王开文，张卫华，等.轮轨接触关系计算方法［J］.交通运输工程学报，2006，6（4）：10-13.

[8] NI Pingtao， WANG Kaiwen， ZHANG Weihua， et al. Calculation Method of Wheel-Rail Contact Relation ［J］. Journal of Traffic and Transportation Engineering， 2006， 6 （4）： 10-13. in Chinese

[9] 倪平涛，刘德刚，曲文强.轮轨与轮轮接触几何计算研究［J］.铁道机车车辆，2012，32（5）：5-9.

[10] NI Pingtao， LIU Degang， QU Wenqiang. Research on Contact Geometry Calculation of Wheel-Rail and Wheel-Roller ［J］. Railway Locomotive & Car， 2012， 32 （5）： 5-9. in Chinese

[11] 倪平涛.轮轨接触几何关系计算再研究［J］.铁道机车车辆，2013，33（4）：7-11，42.

[12] NI Pingtao. Further Research on Wheel-Rail Contact Geometry Calculation ［J］. Railway Locomotive & Car， 2013， 33 （4）： 7-11， 42. in Chinese

[13] 高浩，戴焕云，倪平涛.考虑轮对弹性的轮轨接触点算法［J］.铁道学报，2012，34（5）：26-31.

[14] GAO Hao， DAI Huanyun， NI Pingtao. Algorithm of Wheel-Rail Contact Point for Flexible Wheelset ［J］. Journal of the China Railway Society， 2012， 34 （5）： 26-31. in Chinese

[15] 干锋.高速列车轮轨接触关系研究［D］.成都：西南交通大学，2015.

[16] GAN Feng. Study on Wheel-Rail Contact Relationship of High-Speed Train ［D］. Chengdu： Southwest Jiaotong University， 2015. in Chinese

[17] 翟婉明.车辆-轨道耦合动力学［M］.4版.北京：科学出版社，2015：77-82.

[18] ZHAI Wanming. Vehicle-Track Coupling Dynamics ［M］. 4th ed. Beijing： Science Press， 2015： 77-82. in Chinese

[19] KNOTHE K， GRASSIE S L. Modelling of Railway Track and Vehicle/Track Interaction at High Frequencies ［J］. Vehicle System Dynamics， 1993， 22 （3/4）： 209-262.

[20] 朱丹阳.车桥耦合系统动力相互作用与多点地震响应数值方法研究［D］.大连：大连理工大学， 2015.

[21] ZHU Danyang. Study on Numerical Method of Dynamic Interaction and Multi-Point Seismic Response of Vehicle-Bridge Coupling System ［D］. Dalian： Dalian University of Technology， 2015. in Chinese

[22] 雷晓燕.高速铁路轨道动力学——模型、算法与应用［M］.北京：科学出版社，2015：156-160.

[23] LEI Xiaoyan. High Speed Railway Track Dynamics： Model， Algorithm and Application ［M］. Beijing： Science Press， 2015： 156-172. in Chinese

[24] 吴神花，雷晓燕.交叉迭代算法求解车辆-轨道非线性耦合方程的收敛性讨论［J］.华东交通大学学报，2015，32（3）：23-31.

[25] WU Shenhua， LEI Xiaoyan. Convergence Condition of the Cross-Iterative Algorithm for Vehicle-Track Nonlinear Coupling Equations ［J］. Journal of East China Jiaotong University， 2015， 32 （3）： 23-31. in Chinese

[26] 张斌，雷晓燕，罗雁云.基于Newmark格式的车辆-轨道耦合迭代过程的改进算法［J］.中南大学学报（自然科学版），2016，47（1）：298-306.

[27] ZHANG Bin， LEI Xiaoyan， LUO Yanyun. Improved Algorithm of Iterative Process for Vehicle-Track Coupled System Based on the Newmark Formulation ［J］. Journal of Central South University （Science and Technology）， 2016， 47 （1）： 298-306. in Chinese

[28] 高芒芒，李国龙，杨飞，等.基于强迫振动的列车-轨道-轨下结构垂向耦合动力分析方法及工程应用［J］.中国铁道科学，2021，42（2）：50-58.

[29] GAO Mangmang， LI Guolong， YANG Fei， et al. Vertical Coupling Dynamic Analysis Method and Engineering Application of Vehicle-Track-Substructure Based on Forced Vibration ［J］. China Railway Science， 2021， 42 （2）： 50-58. in Chinese

[30] 雷晓燕，王海.车辆-轨道空间非线性耦合系统交叉迭代算法及应用［J］.振动与冲击，2023，42（10）：136-143.

[31] LEI Xiaoyan， WANG Hai. Cross-Iterative Algorithm and Its Application to the Analysis of Vehicle-Track Spatial Nonlinear Coupling Systems ［J］. Journal of Vibration and Shock， 2023， 42 （10）： 136-143. in Chinese

[32] LEI Xiaoyan， WANG Hai. Dynamic Analysis of the High Speed Train-Track Spatial Nonlinear Coupling System under Track Irregularity Excitation ［J］. International Journal of Structural Stability and Dynamics， 2023， 14 （23）： 1-32.

[33] 张斌.固定辙叉道岔振动系统数值仿真理论及其应用研究［D］.上海：同济大学，2018.

[34] ZHANG Bin. Research on Numerical Simulation Theory of Rigid Frog Turnout and Its Application ［D］. Shanghai： Tongji University， 2018. in Chinese

[35] 侯茂锐.基于多因素耦合作用的高速铁路钢轨廓形偏差限值研究［J］.中国铁道科学，2023，44（6）：125-135.

[36] HOU Maorui. Research on Limit Value of Rail Profile Deviation of High-Speed Railway Based on Multi-Factor Coupling Effect ［J］. China Railway Science， 2023， 44 （6）： 125-135. in Chinese

[37] 卜庆萌，姚林泉.轮轨接触几何关系探讨［J］.苏州大学学报（自然科学版），2011，27（3）：79-84.

[38] BU Qingmeng， YAO Linquan. Study on the Geometric Relation of Wheel-Rail Contact ［J］. Journal of Soochow University （Natural Science Edition）， 2011， 27 （3）： 79-84. in Chinese

[39] 程力.基于车辆-轨道结构垂向耦合系统的数值积分方法的应用与研究［D］.兰州：兰州交通大学，2015.

[40] CHENG Li. Application and Research of Numerical Integration Method Based on Vertical Coupling System of Vehicle-Track Structure ［D］. Lanzhou： Lanzhou Jiaotong University， 2015. in Chinese

[41] 黄龙文.轮轨滚动接触力学分析及疲劳寿命预测方法研究［D］.上海：华东理工大学，2018.

[42] HUANG Longwen. Mechanical Analysis of Wheel-Rail Rolling Contact and Research on Fatigue Life Prediction Method ［D］. Shanghai： East China University of Science and Technology， 2018. in Chinese

[43] 徐金辉.高速车辆-轨道耦合系统随机振动分析及轨道不平顺评价方法研究［D］.成都：西南交通大学，2016.

[44] XU Jinhui. Random Vibration Analysis of High-Speed Vehicle-Track Coupling System and Research on Evaluation Method of Track Irregularity ［D］. Chengdu： Southwest Jiaotong University， 2016. in Chinese

[45] 干锋，戴焕云，高浩.磨耗车轮踏面精确轮轨接触关系计算方法［J］.交通运输工程学报，2014，14（3）：43-51.

[46] GAN Feng， DAI Huanyun， GAO Hao. Calculation Method of Accurate Wheel-Rail Contact Relationship of Worn Abrasion Wheel Tread ［J］. Journal of Traffic and Transportation Engineering， 2014， 14 （3）： 43-51. in Chinese

[47] 雷晓燕.有限元法［M］.北京：中国铁道出版社，2000：177-179.

[48] LEI Xiaoyan. Finite Element Method ［M］. Beijing： China Railway Publishing House， 2000： 177-179. in Chinese

[49] 赫丹，向俊，曾庆元.一种无砟轨道动力学建模的新方法［J］.中南大学学报（自然科学版），2007，38（6）：1206-1211.

[50] HE Dan， XIANG Jun， ZENG Qingyuan. A New Method for Dynamics Modeling of Ballastless Track ［J］. Journal of Central South University （Science and Technology）， 2007， 38 （6）： 1206-1211. in Chinese

基金资助

国家自然科学基金资助项目(51978264)

国家自然科学基金资助项目(52178424)

江西省技术创新引导类计划项目(20223AEI91004)

江西省自然科学重点基金资助项目(20224ACB204018)

江西省教育厅科学技术研究项目(GJJ2205122)

江西省教育厅科学技术研究项目(GJJ171313)

轨道交通基础设施性能监测与保障国家重点实验室自主课题研究项目(HJGZ2023207)