基于 Newmark格式的车辆−轨道耦合迭代过程的改进算法

张斌，雷晓燕，罗雁云(.同济大学 铁道与城市轨道交通研究院，上海，0804；.华东交通大学 铁路环境振动与噪声教育部工程研究中心，江西 南昌，33003)



基于 Newmark格式的车辆−轨道耦合迭代过程的改进算法

张斌1,2，雷晓燕2，罗雁云1
(1.同济大学 铁道与城市轨道交通研究院，上海，201804；
2.华东交通大学 铁路环境振动与噪声教育部工程研究中心，江西 南昌，330013)

摘要：针对车辆−轨道耦合系统振动方程联立求解过程，考虑车辆和轨道2个子系统模型，提出一种将有限元法和非线性接触理论相结合的交叉迭代数值改进算法。该算法将子系统方程非荷载项矩阵进行修正和求逆的预处理，基于Newmark-β积分格式规则，构造具有较高收敛速度及精度的松弛因子函数和收敛准则函数，利用轮轨相互作用力在车辆系统与轨道系统之间的快速交叉迭代，改进并实现轮轨耦合关系的求解。研究结果表明：提出的算法正确、有效，极大地提高了动力学方程数值计算效率；时间步长对系统数值解的稳定性影响显著，松弛因子的合理选择，可起到加速系统迭代和增强迭代稳定性的作用；该算法在解决大型工程振动问题时更具高效求解的优越性。

关键词：耦合动力学；轮轨相互作用；Newmark-β法；交叉迭代过程；数值稳定性

列车−轨道相互作用关系一直是车辆动力学和轨道结构振动研究中最重要的基础性问题，近几十年来，轮轨间作用力引起的机车车辆、轨下基础、桥梁及大地的振动等研究课题十分活跃[1−5]，繁简程度不一的轮轨动力学理论、模型、算法及程序纷纷出现，有效解决了轨道交通发展中一些关键问题，取得了一定的研究成果[6−9]。在各类既有数值分析方法中，动力学方程时程积分求解是重要环节，用以满足轮轨界面形成的车辆系统与轨道系统动态耦合的平衡关系。其中，可采用在每一时间步长内形成车辆−轨道整体动力方程，利用直接积分法统一完成求解，计算稳定性好，但轮轨接触刚度矩阵随车辆位置变化而变化，其在整体刚度矩阵中元素的位置和量值具有时变性，矩阵修正及求逆运算是制约轮轨耦合问题求解效率的根本因素；也可采用分别建立车辆和轨道系统各自的运动平衡方程，以车轮与钢轨之间位移协调或者轮轨力差值作为收敛条件进行方程求解[10−14]，但两子系统间耦合的有效性、收敛过程的可控性等许多问题仍值得进一步深入探讨。在此，本文作者提出一种基于有限元法与非线性接触理论的改进交叉迭代算法，将车辆和轨道两子系统方程非荷载项矩阵进行预处理，同时简化轮轨界面参数传递与协调适应条件，交叉迭代过程控制参数减少至轮轨作用力一项。通过研究时程积分步长、收敛精度对轨道结构振动响应的影响，论证此算法在准确性、稳定性及编程求解效率方面的优越性，为分析机车车辆与线路相互作用问题提供一种高效数值计算方法。

1　车辆系统动力学方程的建立

车辆系统离散为附有二系悬挂的多刚体系统整车模型，如图1所示，取以线路横向对称的半结构研究，考虑车体和转向架的沉浮运动和点头运动，车轮仅考虑沉浮运动，整车模型共有10个自由度，规定位移和力的方向以竖直向上为正，转角和力矩的方向以反时针为正。

图1　整车附有二系悬挂系统车辆模型Fig.1　Vertical model of vehicle with primary and secondary suspensions

图1中，Mc和 Jc分别为车体质量与转动惯量；Mt和Jt分别为转向架质量与转动惯量；Mwi(i=1,2,3,4)为第i个车轮质量；ks1和cs1分别为车辆一系悬挂刚度与阻尼；ks2和Cs2分别为车辆二系悬挂刚度与阻尼；vc和vti(i=1,2)分别为车体和前后转向架的沉浮运动竖向位移；θc和θti(i=1,2)分别为车体和前后转向架的点头运动角位移；vwi(i=1,2,3,4)为第i个车轮的竖向位移；l1和l2分别为车辆固定轴距之半和车辆定距之半。

定义车辆系统的位移向量为

由Hamilton 原理[15]可得车辆系统振动方程

其中：Mu，Cu和 Ku分别为车辆系统质量、阻尼和刚度矩阵；，和 au分别为车辆系统加速度、速度和位移向量；Qu为车辆系统重力向量；Ful为轨道系统对车辆系统的作用力向量。

式中：g为重力加速度。

式中：Fi(i=1,2,3,4)为轮轨相互作用力，可用Hertz接触公式[16]计算如下

其中：vlci和ηi(i=1,2,3,4)分别为第i个轮轨接触处钢轨位移和轨道不平顺值；G 为接触挠度系数，锥形踏面车轮 G 为 4.57r−0.149×10−8(m/N 2/3)，磨耗形踏面车轮G为3.86r−0.115×10−8(m/N 2/3)；r为车轮半径(m)。

2 轨道系统动力学方程的建立

轨道系统亦取半结构研究，板式无砟轨道单元有限元模型如图2所示，钢轨离散为黏弹性点支承梁单元，预制轨道板和水硬性混凝土支承层离散为连续黏弹性支承梁单元。纵向相邻扣件间距为一个单元，单元长度记为l，图2中u1和u4分别为钢轨的竖向位移；θ1和θ4分别为钢轨的转角；u2和u5分别为预制轨道板的竖向位移；θ2和θ5分别为预制轨道板的转角；u3和u6分别为水硬性混凝土支承层的竖向位移；θ3和θ6分别为水硬性混凝土支承层的转角；ky1，ky2和 ky3分别为无砟轨道扣件、CA 砂浆和防水层的刚度；cy1，cy2和cy3分别为无砟轨道扣件、CA砂浆和防水层的阻尼。

图2板式无砟轨道单元模型Fig.2Finite element model of slab track

定义板式无砟轨道单元的节点位移向量为

板式无砟轨道单元的刚度矩阵可表示为

板式无砟轨道单元的质量矩阵可表示为

板式无砟轨道单元的阻尼矩阵可表示为

无砟轨道单元等效节点荷载向量为

其中：Fi计算见式(5)；l=a+b。当车轮不作用于该轨道单元时，为12×1阶零向量。

根据上述轨道单元模型，利用有限元“对号入座法”组集轨道单元质量、阻尼、刚度矩阵和等效荷载向量，可得轨道系统振动方程为

其中：Ml，Cl和Kl分别为轨道系统质量、阻尼和刚度矩阵；，和al分别为轨道系统加速度、速度和位移向量；Ql为轨道系统等效荷载向量。，，，。

3　耦合方程求解算法的改进

研究和实践经验表明[12−14]：采用交叉迭代算法求解车辆−轨道非线性耦合系统振动方程的关键在于迭代过程的稳定性，即在控制精度范围内系统时程响应始终是收敛的。基于Newmark-β积分格式规则的改进交叉迭代法求解耦合系统振动方程的推导过程及步骤如下。

3.1计算准备

按照Newmark积分法格式，对2系统中要求解的二阶常微分方程组进行方程转化。

对于车辆系统运动方程(2)，在t+Δt时刻可写成如下形式：

对轨道系统运动方程(11)，在 t+Δt 时刻可写成如下形式：

其中：积分常数α=0.25，δ=0.5，c0=1/(αΔt 2)，c1=δ/(αΔt)，c2=1/(αΔt)，c3=1/(2α)−1，c4=δ/α−1，c5=Δt(δ/α−2)/2，C6=Δt(1−δ)，c7=δΔt。

将式(12)和式(13)中方程左边部分形成各自等效刚度矩阵，即

式(14)和式(15)在系统求解中具有非时变性，在此完成矩阵求逆，整个数值求解过程中只需要这样1次求逆，将其存储为 Ru和 Rl矩阵，每一时间步直接调用即可。

3.2时间步长循环

首先，从t时刻进入t+Δt时刻时，系统所有变量初值取t时刻的最终值，

令

同时结合式(16)和式(17)，则车辆系统方程式(12)和轨道系统方程式(13)可简化为：

根据积分格式，2 系统相应的速度和加速度可写成如下形式：

从式(20)和式(21)可以看出：车 辆系统和轨道系统之间协调条件的控制参数仅为轮轨作用力一项，很大程度上简化了轮轨界面参数传递过程，2 系统只需进行轮轨相互作用力的交叉迭代求解，即可满足该时刻系统的收敛条件。这里，设在t+Δt时刻，已进行k次迭代，现考察第k+1次迭代：

步骤2：利用松弛法，令

其中：μ为松弛因子，一般取0＜μ＜1，例如μ=0.2～0.4可取得较好效果。

步骤 7：利用轨道位移进行收敛性判别，计算差值：

步骤 8：若收敛性得到满足，则转入下一时刻进行计算，直至整个时域T。

或若收敛性不满足，令 k=k+1，则转入步骤1，继续迭代计算。

从计算步骤可以看出：算法思路清晰易懂，积分过程化繁为简，编程计算效率高。与整体动力方程统一求解方法相比，交叉迭代法避免了每一时间步整体刚度矩阵的重新生成和求逆，尤其当矩阵维数较大时，显著地减少了计算工作量。同时，改进后的算法对振动系统收敛过程的控制效果改善明显，是一种求解车辆−轨道动力学问题的快速数值分析方法。耦合方程求解算法的程序框图如图3所示。，

图3　交叉迭代改进算法计算流程Fig.3　Calculation flow of the improved staggered iterative algorithm

4 算例分析

4.1算法验证

为验证本文提出的交叉迭代改进算法的正确性，与向俊等[17]提出的横向有限条与无砟轨道板段单元的车轨系统竖向振动分析法进行算例比较，计算条件为高速列车(1动+4拖)以200 km/h 速度在板式轨道上运行，取波长为12.5 m、波幅为3 mm的周期性正弦函数为轨道高低不平顺激振源，比较 2种方法的计算结果。采用本文算法计算结果分别如图4～9所示。

图4动车车体垂向加速度时程Fig.4Time history of vertical acceleration of motor vehicle body

图5　动车轮轨垂向力时程Fig.5　Time history of wheel/rail force of motor vehicle

图6　钢轨垂向位移时程Fig.6　Time history of vertical displacement of rail

图7　钢轨垂向加速度时程Fig.7　Time history of vertical acceleration of rail

图8　轨道板垂向位移时程Fig.8　Time history of vertical displacement of slab

图9　轨道板垂向加速度时程Fig.9　Time history of vertical acceleration of slab

从图4～9可以看出：各响应波形符合物理概念，与向俊等[17]计算得到的系统动力响应幅值与变化规律基本一致，证明了本文算法的有效性和可行性。

此外，系统某一时刻(t时刻)拖车轮轨相互作用力和相应轮轨接触处钢轨位移的迭代收敛过程，分别如图10 和图11所示。从图10和图11可见：耦合系统以上一时刻(t−Δt 时刻)平衡状态为初值，经过若干次迭代，逐步在该时刻(t时刻)达到新的平衡状态。

图10　t时刻拖车轮轨力迭代过程Fig.10　Iterative process of wheel/rail force of trailer at t moment

图11　t时刻轮轨接触点钢轨位移迭代过程Fig.11　Iterative process of rail displacement of wheel/railContact point at t moment

4.2时间步长的影响

时间步长Δt分别选取0.05，0.10，0 .50及1.00 ms，考察时间步长对系统方程求解过程的影响。其中：车辆选用CRH3型动车，CRTSⅡ型板式轨道，行车速度为250 km/h，收敛精度设置为1.0×10−7。这里约定，每一时间步内轨道系统方程和车辆系统方程各自分别完成1次求解记为耦合系统1次交叉迭代，每一工况迭代次数按照耦合系统求解过程中最大迭代次数统计。不同时间步长时系统的迭代次数计算结果如表1所示。

表1　不同时间步长系统迭代次数Table1　Number of iterations at different time steps

从表1可以看出：时间步长越小，系统下一时刻的状态越接近上一时刻系统的平衡状态，很容易再次达到平衡位置，因而系统迭代次数少，收敛速度快，且松弛因子对系统求解的影响逐渐降低，但是计算效率低下，此时，配合松弛因子，适当增大时间步长，可有效提高求解效率。例如，计算时间取1s，比较μ=0.3，时间步长Δt从0.5 ms减小到0.1ms，耦合系统求解效率。当 Δt=0.5 ms 时，系统总迭代次数为1/0.000 5×(2×78%+3×22%)=4 440次，其中：2和3为该工况包含的迭代次数，78%和 22%为各自比例，同理，Δt=0.1ms 时，系统总迭代次数为1/0.0001× 2×100%=20 000次，后者为前者的4.5倍，求解效率大大降低，因此松弛因子的使用可以选取更大的时间步长，增强求解时间的经济性。

然而，时间步长取得较大时，系统将严重偏离上一时刻平衡状态，丧失稳定性，导致系统数值解的发散而最终无法收敛。因此，对于相互作用的复杂非线性车轨耦合问题，Newmark数值积分法不再是无条件稳定的，也会产生数值不稳定性现象。

4.3收敛精度的影响

收敛精度ε分别选取1.0×10−7，1.0×10−8，1.0×10−9及1.0×10−10，考虑有砟轨道和无砟轨道2种轨道类型，探讨收敛精度对不同轨道结构动力响应的影响。其中：车辆选用CRH3型动车，轨道参数分别采用铁路干线轨道设计参数[14]和CRTSⅡ型板式轨道参数[6]，行车速度取250 km/h。不同收敛精度对振动系统的影响如图12所示。

从图12可见：无砟轨道耦合系统的迭代次数普遍比有砟轨道的多，原因是无砟轨道单元模型较有砟轨道单元模型位移自由度更多，矩阵维数随之增加引起的，但是在松弛因子位于 0.2～0.4 区间时，两者迭代次数差别不大，且均有明显减少趋势。同时，收敛精度要求越高，系统迭代次数也越多，然而在松弛因子0.3附近，迭代次数并没有随收敛精度提高而增加。因而，松弛因子的合理选择，有助于耦合系统的加速迭代。

此外，收敛精度对系统的动力响应具有显著影响，图13 和图14 所示分别为当收敛精度为1.0×10−7和1.0×10−10时，有砟轨道系统的轮轨相互作用力时程，此时松弛因子取0.1。从图13和图14可见：轮轨力理论值为 70 kN，计算机程序固有相对误差最大为2%～3%，因此，一般动力学分析，收敛精度在1.0×10−7～1.0×10−10之间均可满足计算要求，对于精细化分析则需采用较高的收敛精度。当此工况松弛因子取0.3时，不同收敛精度间系统动力响应差别不明显，反映出松弛因子起增强迭代稳定性的作用。

图12收敛精度对振动系统的影响Fig.12Influence ofConvergence precision on vibration systems

图13　收敛精度1.0×10−7时轮轨作用力时程Fig.13　Time history of wheel/rail force at1.0×10 −7Convergence precision

图14收敛精度1.0×10−10时轮轨作用力时程Fig.14Time history of wheel/rail force at1.0×10 −10Convergence precision

5 结论

1)提出基于 Newmark-β 积分格式的交叉迭代改进数值计算方法。算例对比分析高速列车通过板式轨道时车辆系统和轨道系统的动力响应，验证了算法的正确性。

2)交叉迭代改进算法对耦合方程的积分初值问题和收敛问题引入松弛因子函数，加速系统迭代的同时增强迭代稳定性，较好地反映轮轨相互作用引起的系统动力学响应特征，显著地提高积分求解速度，克服了迭代过程易发散等弊端。

3)迭代收敛精度可根据精细化动力学研究或粗略化工程实际的需要进行设置，该算法在准确性、稳定性和程序求解效率等方面对解决大型工程振动问题具有较强的适用性。

参考文献：

[1]RAY WC,JOSEPH P.Dynamic of structures[M].Berkeley:Computers & Structures,Inc.,1995:120−124.

[2]REZAIEE-PAJAND M,ALAMATIAN J.Numerical time integration for dynamic analysis using a new higher order predictor-corrector method[J].EngineeringComputations,2008,25(6): 541−568.

[3]CHENC,RICLES J M.Stability analysis of direct integration algorithms applied to MDOF nonlinear structural dynamics[J].Journal of Engineering Mechanics,2010,136(4): 485−495.

[4]FERIANI A,MULAS M G,ALIPRANDIC.Time domain iterative procedures for vehicle–bridge dynamic interaction[C]//Proceedings of ISMA2006: InternationalConference on Noise and Vibration Engineering.Belgium: Katholieke Universiteit Leuven,2006:1179−1193.

[5]翟婉明.车辆−轨道耦合动力学[M].3版.北京: 科学出版社,2007:117−124.ZHAI Wanming.Vehicle-trackCoupling dynamics[M].3rd ed.Beijing: Science Press,2007:117−124.

[6]张斌,雷晓燕.基于车辆−轨道单元的无砟轨道动力特性有限元分析[J].铁道学报,2011,33(7): 78−85.ZHANG Bin,LEI Xiaoyan.Analysis on dynamic behavior of ballastless track based on vehicle and track elements with finite element method[J].Journal of theChina Railway Society,2011,33(7): 78−85.

[7]朱志辉,余志武,蒋丽忠,等.高速铁路桥梁及场地土交通振动分析[J].振动工程学报,2012,25(5): 548−555.ZHU Zhihui,YU Zhiwu,JIANG Lizhong,et al.Analysis of bridge-ground vibrations induced by moving loads of high-speed train[J].Journal of Vibration Engineering,2012,25(5): 548−555.

[8]史吏,蔡袁强,徐长节.轨道结构动力响应 Newmark 方法时间积分步长的确定[J].浙江大学学报(工学版),2014,48(2): 228−235.SHI Li,CAI Yuanqiang,XUChangjie.Time step length determination of Newmark method for dynamic responses of railway tracks[J].Journal of Zhejiang University(Engineering Science),2014,48(2): 228−235.

[9]徐庆元,李斌,周小林.高速列车作用下路基上板式无砟轨道动力系数[J].中南大学学报(自然科学版),2011,42(9): 2831−2836.XU Qingyuan,LI Bin,ZHOU Xiaolin.DynamicCoefficient of slab track system on subgrade under high-speed trains[J].Journal ofCentral South University(Science and Technology),2011,42(9): 2831−2836.

[10]张健,谭述君,吴昌华.车辆−轨道非线性耦合动力学的精细积分法及其应用[J].振动与冲击,2012,31(8): 5−10.ZHANG Jian,TAN Shujun,WUChanghua.A precise integration method for vehicle-track nonlinearCoupling dynamics and its application[J].Journal of Vibration and Shock,2012,31(8): 5−10.

[11]张楠,夏禾.基于全过程迭代的车桥耦合动力系统分析方法[J].中国铁道科学,2013,34(5): 32−38.ZHANG Nan,XIA He.A vehicle-bridge interaction dynamic system analysis method based on inter-system iteration[J].China Railway Science,2013,34(5): 32−38.

[12]吴定俊,李奇,陈艾荣.车桥耦合振动迭代求解数值稳定性问题[J].力学季刊,2007,28(3): 405−411.WU Dingjun,LI Qi,CHEN Airong.Numerical stability of iteration scheme for solution of vehicle-bridgeCoupling vibration[J].Chinese Quarterly of Mechanics,2007,28(3): 405−411.

[13]YANG Fuheng,FONDER G A.An iterative solution method for dynamic response of bridge-vehicles systems[J].Journal of Earthquake Engineering and Structural Dynamics,1996,25(2):195−215.

[14]雷晓燕.轨道力学与工程新方法[M].北京: 中国铁道出版社,2002: 42−46.LEI Xiaoyan.New methods in railroad track mechanics & technology[M].Beijing:China Railway Publishing House,2002: 42−46.

[15]陈滨.分析动力学[M].2版.北京: 北京大学出版社,2012: 316−320.CHEN Bin.Analytical dynamics[M].2nd ed.Beijing: Peking University Press,2012: 316−320.

[16]金学松,刘启跃.轮轨摩擦学[M].北京: 中国铁道出版社,2004: 51−57.JIN Xuesong,LIU Qiyue.Wheel/rail tribology[M].Beijing:China Railway Publishing House,2004: 51−57.

[17]向俊,赫丹,曾庆元.横向有限条与无砟轨道板段单元的车轨系统竖向振动分析法[J].铁道学报,2007,29(4): 64−69.XIANG Jun,HE Dan,ZENG Qingyuan.Analysis method of vertical vibration of train and ballastless track system with the lateral finite strip and slab segment element[J].Journal of theChina Railway Society,2007,29(4): 64−69.

(编辑 罗金花)

Improved algorithm of iterative process for vehicle-trackCoupled system based on Newmark formulation

ZHANG Bin1,2, LEI Xiaoyan2, LUO Yanyun1
(1.Institute of Railway and Urban Mass Transit,Tongji University,Shanghai 201804,China; 2.Engineering ResearchCenter of Railway Environment Vibration and Noise,Ministry of Education,EastChina Jiaotong University,Nanchang 330013,China)

Abstract:Based on finite element method and nonlinearContact theory,an improved staggered iterative algorithm for the simultaneous solution process of vehicle-track system equations was presented.In the algorithm,the dynamic model was divided into vehicle subsystem and track subsystem.For the non-load matrices in two subsystems wereCarried out some necessary preprocesses such asCorrection and inversion.Based on rules of the Newmark-β integration scheme,the relaxation factor function andConvergenceCriterion function wereConstructed with higherConvergence speed andConvergence precision respectively.To accomplish the solution of the wheel/railCoupled relationship,the fast staggered iteration of wheel/rail interaction force between vehicle subsystem and track subsystem was applied.The results show that the proposed algorithm isCorrect,effective and able to greatly enhance resolution efficiency of the kinetic equationsCompared withConventional integration mode.The time step has tremendous effect on the stability of numerical solution,and the reasonableChoice of relaxation factor plays aCrucial role in accelerating iteration and enhancing the stability of iterative solution.The proposed algorithm is more efficient in large engineering dynamic response.

Key words:Coupling dynamics; wheel/rail interaction; Newmark-β method; staggered iterative process; numerical stability

中图分类号：U213.2

文献标志码：A

文章编号：1672−7207(2016)01−0298−09

DOI:10.11817/j.issn.1672-7207.2016.01.041

收稿日期：2015−01−29；修回日期：2015−03−29

基金项目(Foundation item)：国家自然科学基金资助项目(U1134107,51208198)；江西省自然科学基金资助项目(20142BAB216001)；江西省教育厅科学技术研究项目(GJJ14393)(Projects(U1134107,51208198)supported by the National Nature Science Foundation ofChina; Project(20142BAB216001)supported by the Natural Science Foundation of Jiangxi Province; Project(GJJ14393)supported by the Science and Technology Program of Education Department of Jiangxi Province)

通信作者：雷晓燕，博士，教授，从事铁路环境振动与噪声控制研究；E-mail: xiaoyanlei2013@163.com