引 言
车辆-轨道耦合系统的动力响应求解是研究车辆振动和轨道结构振动的重要基础［1］.一般情况下，车辆-轨道耦合系统的动力自由度数目较多，其动力学微分方程组的求解只能采用直接数值积分方法.显式数值积分法和隐式数值积分法是直接数值积分法的两种不同形式.其中，隐式数值积分的经典方法有Newmark-β法［2］和Wilson-θ法等，其具有数值积分精度高、稳定性较好的特点，且有利于时间积分步长的选取.
国内外学者对车辆-轨道耦合系统响应的求解算法进行了大量研究［3-10］，提出了多种改进算法及程序，取得了一系列的研究成果.在车辆-轨道耦合系统的数值模拟模型中，轨道结构的长度需要根据模拟时间来选取，模拟时间越长，轨道结构也越长.目前，为了克服这一问题，常采用定点激励，即车辆相对于轨道结构只在一个固定点处振动，同时使轨道不平顺作与列车运行方向相反的相对运动，因此可以极大地缩短轨道结构长度，从而加快计算速率.然而，该算法虽然可以提高计算效率，但并不能模拟列车的真实前进运动状态.文献［2］通过对车辆-轨道垂向耦合系统动力响应计算的迭代过程进行改进，可以实现在一定程度上加快程序运行速度.鉴于此，本文在Newmark-β交叉迭代积分法的基础上，对求解过程进行了改进.本文提出的改进加快算法中的关键步骤在于计算过程中需要不断地对计算模型的车辆位置进行调整，同时，对轨道结构的响应矩阵进行处理，从而实现不增加轨道结构长度而可以增长模拟时间的目的.该算法可以有效地减少轨道结构的轨道单元数量，进而使轨道结构的质量矩阵、阻尼矩阵和刚度矩阵的维数显著降低，为车辆-轨道垂向耦合系统响应求解提供了一种高效精确的算法.
1 车辆-轨道耦合系统垂向振动模型
在车辆-轨道耦合系统垂向振动模型［9］中，一般将整个耦合系统拆分为车辆子系统和轨道子系统，并通过轮对和轨道之间的相互作用力来建立两个子系统的联系.
1.1 车辆子系统
车辆子系统离散为附有二系悬挂的多刚体整车模型，具有10个自由度，分别为车体的垂向运动和点头运动、前后转向架的垂向运动和点头运动以及4个车轮的垂向运动，具体情况如图1所示.
图1 车辆子系统模型
Fig.1 The train subsystem model
图2 轨道子系统模型
Fig.2 The track subsystem model
在图1中，M和J分别为质量和转动惯量，下标c、t和w分别表示车体、转向架和车轮;K和C分别为刚度和阻尼，下标s1和s2分别表示一系悬挂和二系悬挂;z和θ分别为沉浮运动垂向位移和点头运动角位移，下标t i(i=1，2)表示前后转向架，w i(i=1，2，3，4)表示第i个车轮.
1.2 有砟轨道子系统
由于在车辆-轨道耦合系统垂向振动模型中，假定了整个系统模型沿线路方向左右对称.因此，轨道子系统可取轨道对称轴的一半结构来建模，如图2所示.
在图2中，轨道结构离散为四层结构，分别为钢轨、轨枕、道砟和路基.其中，钢轨离散为黏弹性点支承梁单元，轨枕质量和道砟质量均简化为集中质量，K和C分别为刚度和阻尼，下标y1、y2和y3分别表示轨下垫板和扣件、枕下道床的支撑与道砟下路基的支撑.在轨道子系统模型中，一般取两相邻轨枕间的一跨结构作为一个轨道单元，记轨道单元长度为l.
1.3 耦合系统模型的动力学方程
对于车辆子系统和轨道子系统，根据Hamilton原理，得到两子系统的垂向振动方程，其中车辆子系统的垂向振动方程为
式中，M u，C u和K u分别为车辆子系统的质量、阻尼和刚度矩阵;a u，v u和z u分别为车辆子系统的加速度、速度和位移向量;Q u为车辆子系统的重力向量;F u为轨道子系统与车辆子系统之间的相互作用力向量.
有砟轨道子系统的垂向振动方程为
式中，M r，C r和K r分别为轨道子系统的质量、阻尼和刚度矩阵;a r，v r和z r分别为轨道子系统的加速度、速度和位移向量;Q r为轨道结构的重力向量;F r=－F u.
在式(1)和(2)中，车辆子系统和轨道子系统通过轮轨相互作用力向量F u建立联系，其向量表达式为
式中，Fi(i=1，2，3，4)分别为第i个轮轨相互作用力，可用Hertz非线性接触公式来计算.其表达式为［2］
式中，G 为接触挠度系数，本文取 G=3.86r－0.115 × 10－8 m/N2/3，r为车轮半径，m;z w i为第 i个轮对的垂向位移;z rc i和ηi(i=1，2，3，4)分别为第i个轮轨接触点处的钢轨垂向位移和轨道不平顺值，其中位移均以垂直向上为正.
2 耦合方程求解算法的改进
对于车辆-轨道垂向耦合系统模型，若模拟的车辆运行时间越长，则所需的轨道结构也越长，这将导致轨道子系统的质量矩阵、阻尼矩阵和刚度矩阵的维数增大，模拟需要的计算时间将会大幅度增加.为此，本文在Newmark-β交叉迭代数值积分法［11］的基础上提出了一种改进加快算法，相比定点激励法，本文方法可以较为真实地模拟列车在轨道上的前进运行状态，并能求解得到列车位于轨道不同位置时的轨道结构动力响应.
应用Newmark-β交叉迭代数值积分法时，若已知两子系统在t时刻的解为z t，v t，a t，则系统在t+Δt时刻的动力响应可写成如下形式:
式中
其中积分常数 α =0.25，δ=0.5.
为了详细说明改进加快算法的步骤，本文选取了130个轨道单元组成的轨道结构来进行阐述.假设在t时刻已进行了k次迭代，现考察第k+1次迭代.
步骤1 根据第k次迭代结束时的轮对位移
和轮轨接触处的钢轨位移
，根据 Hertz非线性接触公式(4)计算车辆子系统与轨道子系统之间的接触力
步骤2 利用松弛法对轮轨接触力进行修正，令
式中，μ为松弛因子，一般取μ=0.2～0.5可得到较好的效果.
步骤3 将
转换为轮轨接触力向量
并作用于轨道子系统，求解动力方程(2)，得到轨道子系统的位移响应
、速度响应
和加速度响应
.
步骤4 利用步骤3得到的
可得第i个轮轨接触处的钢轨位移
，并由 Hertz非线性接触公式(4)计算轨道子系统对车辆子系统的作用力
.
步骤5 将
作为外荷载作用于车辆子系统，对车辆子系统求解方程(1)，得到车辆系统的位移响应
、速度响应
和加速度响应
.
步骤6 利用步骤5得到的
，计算第i个车轮的位移k
.
步骤7 计算轨道位移差值
式中，
分别为当前迭代步和上一迭代步结束时轨道子系统的节点位移向量.
定义收敛准则:
式中
一般收敛精度 ε 取1.0×10－10～1.0×10－7即可获得较好的模拟结果，本文取ε =1.0 × 10－7.
步骤8 进行收敛性判断，若收敛性不满足，则令k=k+1，进而转入步骤1，继续进行迭代计算;若收敛性满足，则进行下一时间步的计算.
步骤9 当进入t+Δt时间步时，首先判断第4号轮对与钢轨接触点的位置，即作用力F4的位置.若第4号轮对仍位于51号轨道单元内，则按照Newmark-β交叉迭代数值积分法继续进行t+Δt时间步的计算.若第4号轮对已经从51号轨道单元进入到52号轨道单元，如图3所示，此时，根据本文改进加快算法，将车辆子系统相对轨道子系统整体向后移动一个轨道单元，相当于车辆子系统始终位于第51号轨道单元，并使x1=x2，如图4所示.
图3 移动前的轮对位置
Fig.3 The wheelset position before the move
图4 移动后的轮对位置
Fig.4 The wheelset position after the move
对于车辆-轨道垂向耦合系统模型的130个轨道单元，假设t时间步第4号轮对位于第51号轨道单元，此时，轨道单元的位移响应
、速度响应
和加速度响应
分别为
然而，当t+Δt时间步第4号轮对进入第52号轨道单元时，根据Newmark-β交叉迭代数值积分法的原理，t+Δt时间歩第1次迭代的初始值需要使用t时间步轨道子系统的位移响应
、速度响应
和加速度响应
.又由于上述t+Δt时间步计算开始时，将车辆子系统从第52号轨道单元位置后移至第51号轨道单元，则t+Δt时间步第1次迭代的轨道子系统初始位移响应
、速度响应
和加速度响应
需进行如下的处理:
可以发现，为了计算t+Δt时间步的响应，需要使用t时间步的轨道结构位移响应
、速度响应
和加速度响应
，又考虑到t+Δt时间步调整了车辆单元的位置，因此，需要将t时间步的轨道结构响应进行矩阵处理，变换为式(12)～(14)所示的轨道结构位移响应
、速度响应
和加速度响应
.这样，通过轨道结构响应矩阵前移一个单元的方法，即去掉第1个轨道单元响应，尾部增加一个轨道单元响应，保证了车辆子系统后移一个轨道单元的同时，轨道子系统的响应与车辆子系统的位置相对应.此时，即可按步骤1进行t+Δt时间步的计算，如此循环计算下去，直至计算完成.
本文提出的改进加快算法的核心思想，就是在计算过程中不断对车辆子系统的位置和轨道子系统的响应矩阵进行处理，以达到不增加轨道单元数量而可以增长模拟时间的目的.本文算法除了步骤9以外，其他步骤均与Newmark-β算法相同，思路清晰易懂，计算过程简单，程序编写较易，只需要在原先Newmark-β算法的基础上进行少量修改即可完成，而编程计算效率可得到很大的提高.本文提出的改进加快算法通过减少轨道子系统的轨道单元数量，使轨道子系统的质量矩阵、阻尼矩阵和刚度矩阵的维数得到显著降低，为车辆-轨道垂向耦合系统响应求解提供了一种高效的计算方法.
在本质上，文献［2］的算法是对每一个计算时间步的迭代过程进行改进，而本文算法则是对t时间步到t+Δt时间步的传递过程进行改进.这样既可以真实模拟列车在轨道结构上的前进运行状态，又可以不增加轨道结构的长度，只需选取较短的轨道结构，从而能够极大地提高计算效率.
3 算例分析
为了验证本文提出的改进加快算法的精确性和高效性，现与文献［2］中的改进迭代过程算法(传统算法)进行对比分析.轨道不平顺作为车辆-轨道垂向耦合系统的外部激励源，本文采用谱表示-随机函数方法［12-14］模拟的轨道不平顺作为外部激励源，其代表性样本如图5所示.
在算例分析中，车辆子系统选取CRH3型高速列车，列车运行速度为200 km/h，轨道子系统选取60 kg/m有砟轨道结构，其他参数取值见文献［15］.轨道子系统结构由轨道单元组成，其中以两轨枕之间的一跨钢轨、轨枕、道砟和路基作为一个轨道单元，取轨道单元长度l=0.568 m.轨道子系统前后分别取50个轨道单元用以消除边界效应的影响，则轨道结构长度为
图5 轨道不平顺的代表性样本
Fig.5 The representative sample of track irregularity
式中，l2为前后转向架悬挂距离的一半.
因此，整个轨道结构共长73.175 m，需取130个轨道单元.在上述标准参数选取之后，为了考虑模拟时间长度对车辆-轨道耦合系统动力方程求解的影响，考虑以下3种工况:
工况1 模拟时间5 s;
工况2 模拟时间10 s;
工况3 模拟时间15 s.
轨道子系统的轨道单元数量m可由下式计算:
式中，v为车辆运行速度，T为模拟时间，l为轨道单元长度，符号「·表示向上取整数.
图6 工况1下两种算法计算结果比较
Fig.6 Comparisons of the simulation results between the 2 methods in working condition 1
采用MATLAB软件对上述3种工况进行编程计算，表1给出了传统算法和本文算法的计算参数和计算耗时结果比较;图6～8给出了传统算法和本文算法的计算结果比较.
表1 各工况响应计算时间对比
Table 1 Comparison of response calculation time of all working conditions
working condition working condition 1 working condition 2 working condition 3 number of elements in the traditional algorithm N T 619 1 108 1 597 number of elements in the proposed algorithm N P 130 130 130 calculation time in the traditional algorithm t T/s 2 819 16 210 43 660 calculation time in the proposed algorithm t P/s 298 520 753 calculation time ratio S=t T/t P 9.46 31.17 57.98
图7 工况2下两种算法计算结果比较
Fig.7 Comparisons of the simulation results between the 2 methods in working condition 2
图8 工况3下两种算法计算结果比较
Fig.8 Comparisons of the simulation results between the 2 methods in working condition 3
从图6～8可以看出，本文算法与传统算法的计算结果几乎完全一致，表明本文算法具有良好的计算精度.同时，从表1中可知，随着模拟时间的增加，传统算法的轨道子系统质量矩阵、阻尼矩阵和刚度矩阵的大小都呈现指数倍增长，使得计算时间也呈现出指数倍增长.而本文算法由于单元数量保持不变，因而计算时间仅呈小幅度增加，从而极大地提高了本文算法的计算效率.
4 结 论
以Newmark-β交叉迭代数值积分法为基础，本文提出了一种求解车辆-轨道垂向耦合系统动力响应的改进加快算法.在计算过程中不断地定位车辆子系统的位置，通过调整车辆子系统的位置和轨道子系统的动力响应矩阵，可以在保证真实模拟车辆在轨道结构上前进运行状态的情况下，不额外增加轨道子系统的单元数量.研究表明，本文提出的改进加快算法只需要选取较短的轨道结构长度，即可对不同动力响应持时的工况进行模拟，既能保证模拟结果的精确性，又能极大地提高计算效率，为求解车辆-轨道垂向耦合系统的动力响应提供了一种高效精确的算法.
［1］ 徐磊，翟婉明.轨道结构随机场模型与车辆-轨道耦合随机动力分析［J］.应用数学和力学，2017，38(1):67-74.(XU Lei，ZHAI Wanming.The random field model for track structures and vehicle-track coupled stochastic dynamic analysis［J］.Applied Mathematics and Mechanics，2017，38(1):67-74.(in Chinese))
［2］ 张斌，罗雁云，雷晓燕.改进迭代过程的车轨耦合振动数值解法及应用［J］.工程力学，2016，33(3):128-134.(ZHANG Bin，LUO Yanyun，LEI Xiaoyan.Numerical approach for vehicle-track coupling vibration based on improved iterative process and its application［J］.Engineering Mechanics，2016，33(3):128-134.(in Chinese))
［3］ YANG F，FONDER G A.An iterative solution method for dynamic response of bridge-vehicles systems［J］.Earthquake Engineering ＆ Structural Dynamics，1996，25(2):195-215.
［4］ 娄平，曾庆元.移动荷载作用下板式轨道的有限元分析［J］.交通运输工程学报，2004，4(1):29-33.(LOU Ping，ZENG Qingyuan.Finite element analysis of slab track subjected to moving load［J］.Journal of Traffic and Transportation Engineering，2004，4(1):29-33.(in Chinese))
［5］ JO J S，JUNG H J，KIM H J.Finite element analysis of vehicle-bridge interaction by an iterative method［J］.Structural Engineering and Mechanics，2008，30(2):165-176.
［6］ 雷晓燕，张斌，刘庆杰.列车-轨道系统竖向动力分析的车辆轨道单元模型［J］.振动与冲击，2010，29(3):168-173.(LEI Xiaoyan，ZHANG Bin，LIU Qingjie.Model of vehicle and track elements for vertical dynamic analysis of vehicle-track system［J］.Journal of Vibration and Shock，2010，29(3):168-173.(in Chinese))
［7］ ZHANG J，GAO Q，TAN S J，et al.A precise integration method for solving coupled vehicletrack dynamics with nonlinear wheel-rail contact［J］.Journal of Sound and Vibration，2012，331(21):4763-4773.
［8］ SUN W，ZHOU J，THOMPSON D，et al.Vertical random vibration analysis of vehicle-track coupled system using Green’s function method［J］.Vehicle System Dynamics，2014，52(3):362-389.
［9］ 吴神花，雷晓燕.交叉迭代算法求解车辆-轨道非线性耦合方程的收敛性讨论［J］.华东交通大学学报，2015，32(3):23-31.(WU Shenhua，LEI Xiaoyan.Convergence condition of cross iterative algorithm for vehicle-track nonlinear coupling equations［J］.Journal of East China Jiaotong University，2015，32(3):23-31.(in Chinese))
［10］ 孙宇，翟婉明.基于格林函数法的车辆-轨道垂向耦合动力学分析［J］.工程力学，2017，34(3):219-226.(SUN Yu，ZHAI Wanming.Analysis of the vehicle-track vertical coupled dynamics based on the Green’s function method［J］.Engineering Mechanics，2017，34(3):219-226.(in Chinese))
［11］ 翟婉明.车辆-轨道耦合动力学［M］.北京:科学出版社，2015.(ZHAI Wanming.Vehicle-Track Coupling Dynamics［M］.Beijing:Science Press，2015.(in Chinese))
［12］ 刘章军，方兴.平稳地震动过程的随机函数-谱表示模拟［J］.振动与冲击，2013，32(24):6-10.(LIU Zhangjun，FANG Xing.Simulation of stationary ground motion with random functions and spectral representation［J］.Journal of Vibration and Shock，2013，32(24):6-10.(in Chinese))
［13］ LIU Z J，LIU W，PENG Y B.Random function based spectral representation of stationary and non-stationary stochastic processes［J］.Probabilistic Engineering Mechanics，2016，45:115-126.
［14］ 刘章军，张传勇，齐东春.高速列车乘坐舒适度评价的概率密度演化方法［J］.振动与冲击，2018，37(14):149-155.(LIU Zhangjun，ZHANG Chuanyong，QI Dongchun.Probability density evolution method for the ride comfort evaluation of high speed trains［J］.Journal of Vibration and Shock，2018，37(14):149-155.(in Chinese))
［15］ 雷晓燕.高速铁路轨道动力学:模型、算法与应用［M］.北京:科学出版社，2015.(LEI Xiaoyan.High Speed Railway Track Dynamics:Model，Algorithm and Application［M］.Beijing:Science Press，2015.(in Chinese))