基于传递矩阵法的车辆-轨道垂向耦合系统动力分析方法

朱志辉， 王盈莹， 龚 威， 冯乾朔

(1.中南大学 土木工程学院，长沙 410075； 2.高速铁路建造技术国家工程实验室，长沙 410075)

铁路车辆-轨道耦合系统由移动车辆和连续轨道结构组成。列车沿轨道高速运行将引起轨道结构振动，强烈的振动不仅会直接影响轨道的工作状态和使用性能，还会降低列车的行车平稳性和安全性，影响乘客乘坐舒适度[1]。因此，准确和有效地评估车辆和轨道间的相互作用是至关重要的。

建立准确的轨道结构动力学模型是车辆-轨道耦合动力学理论研究的重要部分。有限元直接刚度法(direct stiffness method, DSM)可直接基于轨道结构有限元模型建立其动力方程[2-4]。这种方法可准确计算轨道结构的局部与高频振动[5]，但轨道模型的自由度数目会随着轨道长度的增长而增加，当涉及较长轨道结构动力分析时，不可避免会产生大量结构自由度，使得动力方程系数矩阵规模庞大，导致建模和计算的困难[6]。为了克服这个问题，翟婉明[7]假设车辆相对于轨道结构只在一个固定点振动，提出通过加载车轮上的轨道随机不平顺性得到车辆响应的近似办法，具有较高的计算效率，但忽略了轨道结构变化的影响。Guo等[8-10]采用模态叠加法建立轨道结构动力方程，缩减了轨道自由度数目，但由于轨道结构的局部高频振动特性，导致其主导振型难以选取[11]。移动轨道技术[12-13]通过删除非荷载作用处轨道结构有效降低了自由度数目，但新增轨道连接处存在位移不协调问题，新增单元处会出现动力响应突变[14]。

传递矩阵法(transfer matrix method，TMM)的基本思想是将周期性结构分割为若干元胞结构，将元胞的力学特性用传递矩阵表示，建立元胞结构状态矢量传递模型，实现整个周期性结构的求解。传递矩阵法基于元胞结构进行分析，无需建立结构整体动力方程，避免了大型线性方程组的求解[15-16]，降低了计算过程中的内存占用。Mead等[17-21]采用传递矩阵法分别针对周期性的弹性支承梁、框架结构、薄壁构件、拱桥和钢桥进行了研究，取得了良好的效果。对于链式系统，状态矢量的传递仅相当于元胞单元传递矩阵相乘，传递矩阵连续相乘即可得到系统总传递矩阵。然而，由于单元传递矩阵连续相乘，累积误差被加剧，可能产生数值不稳定现象[22]。针对该问题，Horner等[23]结合Riccati变换与传递矩阵法，提出了一种Riccati传递矩阵法，通过传递矩阵的递推建立系统传递关系，有效提高了传递矩阵法的数值稳定性。Xue[24]在元胞结合面处引入刚度方程假定，将传递矩阵连续相乘演变为元胞结合面刚度方程的递推传递，有效减小了因矩阵连续相乘产生的累积误差。

轨道结构具有明显的周期特性，Ma等[25-26]利用这一特性，有效地解决了轨道结构的频域响应分析问题，但该类方法无法考虑轨道的非线性特征和轨道构成差异性。因此，本文利用传递矩阵方法建立轨道结构传递模型并求解其动力响应，从而实现不增加轨道系统计算自由度数目而可以增长轨道结构计算长度的目的。同时，当轨道结构纵向构成有差异，如出现轨枕空吊、扣件失效等情况时，只需将该处元胞更换为包含差异构成的元胞，不影响整体传递关系的建立。

本文将车辆-轨道耦合系统分解为车辆子系统和轨道子系统，两子系统间通过轮轨相互作用力实现耦合。其中，车辆子系统采用多刚体动力学方法建立动力方程；对于具有显著周期性特征的轨道子系统，以有砟轨道和CRTSⅡ型板式无砟轨道为例，基于TMM将轨道结构的周期性重复部分划分为若干类元胞结构，建立元胞结构的动力方程，并基于引入刚度方程假定的Riccati传递矩阵法，确定元胞内部及相邻元胞间的传递关系，从而建立周期性轨道结构传递矩阵模型。通过对比TMM与DSM计算的CRH2型动车通过无砟轨道结构引起的车辆和轨道结构动力响应结果，验证了本文方法的可行性和高效性。

1 轨道结构传递矩阵模型

根据钢轨下部有无道砟支撑层，现有轨道结构可分为有砟和无砟轨道两类，其结构形式和力学特征均有所不同，因此采用传递矩阵法对其进行动力分析时，应针对不同类型的轨道结构，提出相应的元胞划分方案。CRTSⅡ型板式无砟轨道是我国广泛应用的无砟轨道结构形式之一，京津城际、京沪高铁等均采用了此种轨道结构[27]。故本文以有砟轨道和CRTSⅡ型纵连板式无砟轨道为例，提出元胞划分方案，推导元胞间的动力传递关系，建立基于传递矩阵法的轨道结构元胞模型。轨道系统左右对称，为简便起见，取半结构进行介绍，且只考虑在竖向平面内的动力响应。

1.1 元胞结构的划分

轨道结构由重复的单跨轨道构成，因此应用传递矩阵法时，仅需对单跨轨道进行划分。为便于建立元胞内部传递关系，元胞结构的划分应遵循仅元胞端部存在自由度，且相邻元胞交界面自由度数目相同的原则。

有砟轨道垂向整体模型如图1(a)所示，其中钢轨视作连续弹性离散点支撑的垂向Timoshenko梁模型，考虑竖向位移Zr及转动角θr，长度为L；轨枕(及道砟质量)等效为集中质量ms，考虑竖向位移Zs，相邻两轨枕间距为l；扣件和道砟等效为弹簧-阻尼器，考虑竖向刚度系数kf，kb及阻尼系数cf，cb。为对钢轨进行较为精细的划分，如图1(b)所示，将有砟轨道单跨轨道划分为3类元胞结构：元胞A为左侧支撑结构，由钢轨、轨枕、扣件和道砟组成；元胞B为中间连接结构，由钢轨组成；元胞C为右侧支撑结构，构成与元胞A相同，由于相邻两元胞包含同一支承结构，故二者支承结构参数均为整体的1/2(端部元胞除外)。

CRTSⅡ型板式无砟轨道如图2(a)所示，钢轨、轨道板和底座板被视为连续弹性离散点支撑的Timoshenko梁模型，长度为L，考虑竖向位移Zr、Zt、Zc及转动角θr、θt、θc,下标r、t、c分别表示钢轨、轨道板、底座板；扣件、CA砂浆和基础等效为周期性离散分布的黏滞阻尼和线性弹簧，间距为l，考虑竖向刚度系数kf、kc、ks和阻尼系数cf、cc、cs，下标f、c、s分别表示扣件、CA砂浆、基础。类似于有砟轨道，将CRTSⅡ型板式无砟轨道单跨轨道结构划分为3类元胞结构，如图2(b)所示。其中，元胞A和C均由钢轨、轨道板、底座板、扣件、CA砂浆和地基组成，由于相邻两元胞包含同一支承结构，故二者支承结构参数均为整体模型的1/2(端部元胞除外)；元胞B由钢轨、轨道板和底座板组成。

(a) 整体垂向模型

(b) 单跨轨道划分示意图图1 有砟轨道垂向模型划分示意图Fig.1 Structure division of ballast track

(a) 整体垂向模型

(b) 单跨轨道划分示意图图2 CRTSⅡ型板式无砟轨道垂向模型划分示意图Fig.2 Structure division of CRTSⅡ slab ballastless track

1.2 元胞结构动力方程

任意取整体结构中第i个轨道元胞，其运动方程可表示为

(1)

采用Newmark-β直接积分法[29]对元胞动力方程进行求解，式(1)可写为

(2)

图3 元胞i输入端到输出端传递关系示意图Fig.3 The transmission of vectors from the input to the output of cell i

对于元胞i

(3)

如图2所示，元胞i输入和输出端自由度对应的位移和内力可表示为

(4)

式中:Z、θ、Q和M分别表示位移、转角、剪力和弯矩;下标r、t和c分别表示钢轨、轨道板和底座板。

(5)

1.3 传递模型的建立

由位移和内力的连续条件[31]，元胞结合面处的状态矢量传递关系可表示为

Ui-1,O=Ui,I,Ni-1,O=-Ni,I

(6)

引入类似于状态矢量的广义Riccita变换，假设元胞i左侧结合面处内力向量与位移向量的广义刚度方程为

Ni,I=SiUi,I+Ei(i≥2)

(7)

式中:Si是结合面处刚度方程的系数矩阵;Ei是结合面处的等效外力向量。

展开式(5)，可得：

(8)

(9)

将式(6)和式(7)代入式(8)，可得由第i个元胞输出端位移响应求解第i-1个元胞输出端响应的传递关系式

(10)

将式(10)代入式(9)，可得元胞i+1输入端内力与位移向量的广义刚度方程为

Ni+1,I=Si+1Ui+1,I+Ei+1

(11)

其中：

(12)

(13)

式(12)和式(13)即为相邻剖面间，刚度方程的系数矩阵和等效外荷载向量的传递关系。

式(11)中，令i=2，有：

N2,I=S2U2,I+E2

(14)

将式(6)代入式(14)，则上式可表示为

N1,O=-S2U1,O-E2

(15)

对于元胞1，扩展式(5)可得：

(16)

(17)

轨道结构两端无内力作用，N1,I是已知的，因此由式(16)可得元胞1的输入端位移为

(18)

将式(18)代入式(17)可得元胞1的输出端内力为

(19)

对照式(15)可知

(20)

(21)

将S2和E2代入式(12)和式(13)，可完成Si和Ei在整个结构中由输入端至输出端作刚度方程的连续传递，进而建立整体结构的状态矢量传递模型。考虑结构输出端边界条件，通过式(7)可得结构输出端的位移和内力响应，将其代入输出端响应传递关系式(10)，由输出端至输入端依次求解，可得各元胞的位移响应，并通过Newmark-β方法解得其速度和加速度响应。

由上述推导过程可知，采用传递矩阵法对轨道子系统进行动力分析时，无需建立轨道结构整体模型，仅需拆分其单跨轨道为若干元胞结构，并依据元胞结构间的传递关系，建立轨道结构的整体传递模型即可。对于垂向有砟/无砟轨道模型，无论轨道整体长度如何变化，计算过程中所涉及的最大矩阵阶数始终为3/6。

2 车辆-轨道耦合系统动力分析模型

本文在车辆-轨道动力分析模型中，将整个耦合系统划分为车辆子系统和轨道子系统，通过轮轨间的相互作用力来建立两个子系统的联系。其中，轨道子系统模型通过轨道结构传递矩阵模型建立。

2.1 车辆子系统

车辆视作多刚体模型，单节车辆由1个车体、2个转向架、4个轮对、4个一系悬挂系统和2个二系悬挂系统组成，如图4所示。车体、转向架和轮对采用刚体模拟，悬挂系统通过线性弹簧-阻尼器模拟。本文仅关注车辆在竖向平面内的运动，考虑车体的沉浮zc和点头βc、转向架的沉浮zt1、zt2和点头βt1、βt2，以及4个轮对的沉浮zwj1～zwj4，共10个自由度。图中：mc、mt和mw别为车体、转向架和轮对的质量；kpz、cpz和ksz、csz分别为车辆一系和二系悬挂系统的刚度和阻尼系数；Icy和Ity分别为车体和转向架的点头转动惯量；lc为车辆定距之半；lt为车辆轴距之半；v为车辆运行速度。车辆主要参数参见文献[32]。

车辆子系统的动力方程可表示为

(22)

图4 车辆模型Fig.4 Vehicle model

2.2 轨道子系统

采用TMM建立轨道子系统模型，无需建立轨道结构整体模型，仅需将周期性轨道结构划分为若干元胞结构，建立元胞结构的动力方程，并依据元胞结构间的传递关系，建立轨道结构传递模型进行求解。对于元胞m，其动力方程可表示为

(23)

(24)

式中:Brvm(1×(Nm,INm,O))为布尔矩阵，该矩阵中与元胞结构钢轨节点对应自由度的元素为1，其余元素为0;Nm,I和Nm,O分别为元胞m输入端和输出端的自由度数目;Fr为轨道子系统的外荷载向量。

2.3 轮轨间相互作用力

车辆子系统和轨道子系统之间的动力相互作用通过轮轨接触实现，本文采用基于切线斜率法的线性Hertz接触模型模拟轮轨接触关系。切线斜率法是指过非线性赫兹接触曲线中静态轮轨力P0对应的点作切线，切线的斜率即为轮轨接触弹簧刚度kh，可根据下式求解

(25)

式中：G为轮轨接触常数,m/N2/3；P0为静轮重,N。

在式(22)和(24)中，车辆子系统和轨道子系统通过轨道不平顺激励和轮轨间压缩引起的轮轨力实现耦合。需要注意的是，假设车辆在平面内运动，轨道为三维空间模型。因此，每个轮对由两个相同的赫兹弹簧连接在两条钢轨上。

式(22)中的车辆外荷载向量Fv可写为

Fv=Fvg+Fvi+Fvr

(26)

其中，

Fvgj=[mcg0mtg0mtg0mwgmwgmwgmwg]T

(27)

Fvij=[0 0 0 0 0 0fvij1fvij2fvij3fvij4]T

(28)

fviji=2khr(xji)

(29)

式中:Nv为车辆子系统的列车节数；Fvgj(10×1)为第j节列车所受的重力荷载，由各刚体所受重力组成；Fvij(10×1)为车辆所受的轨道不平顺激励；fviji为第j节车辆第i个轮对受到的轨道不平顺激励；r(x)表示轨道不平顺取值，由轮对在轨道上的位置决定。

Fvr(Nv×1)为车辆所受的由轮轨间压缩引起的轮轨力，需结合轮下钢轨进行求解，如图5所示。

图5 车辆和轨道子系统耦合方式Fig.5 Vehicle and track subsystem coupled mode

对于第j节车的第i个轮对，fvrji可表示为

fvrji=kh(zrwji-zwji)

(30)

式中:(zrwji-zwji)表示第j节车的第i个轮对所在位置的轮轨压缩量;zwji为轮对竖向位移;zrwji为轮下钢轨的竖向位移。可由轮下元胞中钢轨单元节点的位移插值求得

zrwji=XmNT(x)|x=xji

(31)

其中，

Xm=Brvm[Um-1,OUm,O]T

(32)

式中:NT(x)为钢轨单元的形函数向量转置；Xm为轮下钢轨单元的位移；Brvm同上；Um-1,O和Um,O分别为位移矩阵中与轨下元胞的输入端和输出端对应的元素。

式(24)中的轨道外荷载向量Fr为

Fr=Fri+Frv

(33)

其中，

(34)

(35)

式中:Fri为轨道所受的不平顺激励；Frv为轨道所受的由轮轨间压缩引起的轮轨力；Njis(1×Nr)为第j节车第i个轮对所在的第s个钢轨单元(1个轮对与2个钢轨单元存在接触)的形函数。

2.4 计算流程

图6给出了基于本文方法的车辆-轨道耦合系统动力响应计算分析流程图，具体计算步骤为：①采用多刚体动力学方法建立车辆动力方程；②将轨道结构划分为元胞结构，建立元胞有限元模型，并采用TMM建立轨道结构传递矩阵模型；③令系统初始响应为零；④根据列车位置确定轨下元胞，计算轮对和钢轨位移差，并通过轮轨接触关系和轨道不平顺求解式(22)和(24)中轮轨间相互作用力；⑤施加轮轨力至对应的轮对和钢轨自由度，通过数值积分法和传递矩阵法求解当前时间步车辆及轨道的动力响应；⑥判断收敛性；⑦存储响应，并将其作为下一时间步的初始状态，重复步骤④和⑤直至计算结束。

图6 基于TMM的车辆-轨道耦合系统动力响应计算流程图Fig.6 The calculation flow of vehicle-track coupled system dynamic response based on TMM

3 方法验证

为验证上述关于车辆-轨道耦合系统传递矩阵算法的正确性，以高速列车作用下的无砟轨道结构为例，对比了TMM和DSM计算所得的轮轨力、车辆位移和轨道动力响应结果。并通过两种方法计算耗时的对比，进一步验证本文方法的高效性。

CRTSⅡ型板式无砟轨道结构如图7所示，轨道长度为60 m，扣件间距为0.6 m；钢轨通过扣件直接与轨道板(内置轨枕的预制板)相连，对称分布于距轨道板中心0.745 m处；轨道板长6 m，宽2.55 m，高0.2 m，其下铺设有CA砂浆层；底座板长6 m，宽3.25 m，高0.3 m，通过板下的滑动层与下部结构相连。假定单节CRH2动车以300 km/h的速度沿钢轨从左至右移动，计算时考虑我国新规范高速铁路无砟轨道不平顺谱样本，结构阻尼比取为0.02，数值积分步长取1/10 000 s。

图7 CRTSⅡ型板式无砟轨道结构示意图Fig.7 CRTSⅡslab ballastless track

对于空间的轨道结构，其元胞划分方案与平面轨道结构类似，其中元胞A、B和C长度均为0.1 m。采用有限元法建立轨道结构的元胞模型，并导出其质量、刚度和阻尼矩阵进行元胞动力方程求解。元胞有限元模型中，钢轨采用空间梁单元Beam188模拟，轨道板和底座板采用空间板单元Shell181模拟，扣件、砂浆层和滑动层均通过弹簧-阻尼器单元模拟，具体建模参数参考文献[34]。

图8、图9和图10给出了TMM和DSM计算所得的第一轮对轮轨力、第一节车体响应及轨道跨中截面处钢轨响应时程曲线。其中，虚线表示采用TMM的计算结果，实线表示DSM的计算结果，从图中可以看出，两种方法的计算结果有很好的一致性。以DSM计算结果为标准，分析两种方法计算所得各项响应的最大值及相对误差，如表1所示。其中最大值误差均小于1%，表明基于本文方法建立的轨道结构传递矩阵模型，应用于车-轨耦合系统动力响应分析时具有较高的计算精度。

图8 第一轮对轮轨力时程曲线Fig.8 Time history curve of vertical wheel-rail contact force of the first wheel set

图9 第1节车体竖向加速度时程曲线Fig.9 The dynamic response time history curve of the vehicle

(a) 竖向位移时程

(b) 竖向加速度时程图10 跨中截面处钢轨动力响应时程曲线Fig.10 The dynamic response time history curve of midpoint of rail

表1 两种方法计算精度对比Tab.1 Comparison of calculation accuracy

为进一步说明本文方法具有更高的计算效率，对比了TMM和DSM两种方法的计算耗时，如表2所示。其中，采用DSM求解轨道结构动力响应时，需要建立无砟轨道结构的整体模型，动力方程最大矩阵阶数为50 484；采用TMM时，仅需建立元胞结构模型，极大地减少了轨道自由度数目，动力方程最大矩阵阶数仅为84，为整体模型的1/601。对比TMM和DSM的计算耗时，当无砟轨道长度为60 m时，TMM计算耗时为DSM的17.1%。表明基于本文方法求解车辆-轨道耦合系统模型动力响应时，仅需对各个轨道元胞结构的动力方程进行求解，无需建立和求解整体轨道结构的动力方程，有效地降低了计算过程中动力方程系数矩阵的规模，提升了计算效率。

表2 两种方法计算耗时对比

4 结 论

本文针对车辆-轨道耦合系统动力分析问题，利用轨道结构的周期性特征，基于传递矩阵法提出了一种便捷的建模及求解方案。针对有砟轨道和CRTSⅡ型无砟轨道系统结构，分别提出了对应的元胞划分方案和传递矩阵模型，并应用于车辆-轨道耦合系统动力响应分析。以高速动车通过无砟轨道结构为例，对其正确性进行了验证。得到了以下结论：

(1) 本文方法计算结果与直接刚度法计算结果吻合度极高，表明基于本文方法应用于车-轨耦合系统动力响应分析时具有较高的计算精度。

(2) 本文方法在建立车辆-轨道耦合系统模型时应用了轨道结构的周期特性，只需建立其单跨轨道划分而成的元胞模型，无需建立其整体模型，使得建模更为便捷。

(3) 基于传递矩阵法求解车辆-轨道耦合系统模型动力响应时，仅需对各个轨道元胞结构的动力方程进行求解，无需组装整体轨道结构的动力方程，有效地降低了计算过程中涉及矩阵的规模，提升了计算效率。