基于虚拟激励法列车脱轨系数概率统计研究

周子骥，张 楠，孙琪凯

(北京交通大学土木建筑工程学院，北京 100044)

在车-桥耦合系统动力响应研究中，轨道不平顺是引起耦合系统振动的重要激励源。由于轨道不平顺随机函数是平稳Gauss 随机过程，因此引起的车桥耦合振动也是典型的随机过程[1]。确定性分析只针对给定参数得到的单一结果，难以体现车桥系统安全性和平稳性指标的概率分布规律，必然存在一定的风险，由此开启了基于车桥随机振动的桥梁安全性研究。

传统车-桥随机振动研究常采用的Monte Carlo法，通过谐波叠加等方法将轨道不平顺功率谱生成多组不平顺时域样本，根据不同的时域样本计算车-桥系统动力响应，再采用Monte Carlo 法统计计算结果。此方法是基于大样本容量的随机分析方法[2]。Chatterjee 等[3]根据桥面随机不平顺激励样本，采用基于Monte Carlo 法计算桥梁冲击系数的均值和偏差。虽然Monte Carlo 法原理简单，在车-桥随机振动研究中被广泛使用，但是由于其必须取足够数量的样本才能保证统计结果的可靠性，效率极低[4]。在此基础上，学者们又进一步运用重要抽样法、一次二阶矩等方法研究桥梁结构可靠度[5−6]。同时，部分学者通过改进迭代计算方法来提高计算的精确性与时效性[7]。随着有限元仿真的大力发展，学者们利用有限元模拟建立车-轨系统研究了该系统的随机振动特性[8]。为进一步降低随机振动分析的计算量，同时提高随机振动分析的准确性，林家浩等[9]提出了随机振动的虚拟激励法。

基于虚拟激励法，余志武等[2,10]、李小珍等[11−12]、赵岩等[13]、何旭辉等[14]等将其用于车-桥随机系统，研究了列车荷载作用下车-桥耦合系统的随机振动特性，并采用3 σ法则对桥梁响应上下限值进行估计，研究了桥梁的随机动力特性。王立东等[15]提出了一种考虑轨道随机不平顺的列车地面振动预测方法，将土体简化为2.5D 有限元匹配模型，研究了列车随机振动对土体的影响。在此基础上，部分研究亦在车桥系统分析中的同时考虑了其他类型的随机荷载，分析了随机荷载对桥梁动力响应极值的影响[16−17]。但是大多数研究主要针对桥梁结构安全性，对于列车运行安全指标的研究很少。在结构安全方面的研究，常用方法是假设概率分布函数得到安全界限值，然后开展讨论[18]。对于车-桥系统而言，列车安全事故发生概率要高于因列车运行造成桥梁结构破坏概率，而且脱轨系数不服从正态分布，需采用新的概率分布函数进行统计研究。因此对于列车行车安全指标的研究有较大的意义。

基于上述研究，本文采用竖向密贴横向线性蠕滑轮轨关系建立车-桥动力学模型，采用虚拟激励法求得列车轮轨力功率谱；运用谐波叠加法，通过每一时刻三角级数分解求非平稳轮轨力功率谱的时域样本；再由极值理论对列车脱轨系数进行安全统计；最后选取置信度99.73%的脱轨系数极值与多次时域样本最大值进行比较，验证使用极值理论研究脱轨系数的正确性。

1 车-桥耦合系统计算模型与轮轨关系

将车辆和桥梁看成两个子系统，分别建立车辆和桥梁系统的模型。

根据拉格朗日方程可建立车辆系统动力学方程，列车运动方程为：

式中：Mv、Cv和Kv分别为车辆系统的质量矩阵、阻尼矩阵和刚度矩阵；Fv为车辆系统受到的车桥系统相互作用力；X、X˙ 、X¨分别是位移向量、速度向量和加速度向量。车辆单元位移向量的自由度排列顺序为：

式中：下标c 表示车体；t1、t2 依次表示前、后转向架；w1、w2 依次表示与前转向架相连的两轮对；w3、w4 依次表示为与后转向架相连的两轮对。

在车桥耦合系统相互作用研究中，为降低计算工作量，常采用振型叠加法建立桥梁动力方程。桥梁系统模型采用空间梁单元形式，每个单元节点有6 个自由度，分别为沿着x轴、y轴、z轴的平移和绕着x轴、y轴、z轴的转动。本文不考虑轨道结构，假设轨道结构和桥梁没有相对位移，轮轨力通过轨道直接作用于桥梁上。基于此可以得到桥梁运动方程为：

式中：Xb是桥梁广义坐标向量；w是频率向量；Φ为振型矩阵。

本文轮轨相互作用模型采用横向线性蠕滑，垂向密贴的轮轨关系假定。轮轨垂向假定即假定车辆子系统坐标z方向，轮轨间无相对位移，此方向的相互作用力可以由轮轨间该方向的相对运动状态确定；横向线性蠕滑理论即满足虚拟激励法对线性系统的要求，又可考虑轮轨间横向运动和相应作用力之间的线性关系。本课题组研究发现，对于我国常见的车轮半径，线性蠕滑参数与非线性蠕滑参数误差在0.13%～4.43%。因此采用横向线性蠕滑理论可以满足工程精度要求[19]。

2 基于虚拟激励法的脱轨系数统计研究

2.1 车-桥耦合时变系统随机振动分析

由第1 节可知，列车子系统和桥梁子系统通过轮轨接触关系组成整体耦合时变系统，其耦合时变系统动力方程：

式中：角标v 和b 分别为车辆子系统和桥梁子系统；M、K和C分别为质量矩阵、刚度矩阵和阻尼矩阵；F为力向量。

列车桥梁耦合系统所受外荷载F(t)包括车辆自重引起的确定性动荷载Fg和轨道不平顺引起的随机性荷载Fr(t)。即：

尽管车桥子系统都是线性系统，但是由于列车在桥上的位置是随时间不断变化的，因此车桥系统具有非平稳特性。可将随机激励F˜r(t)表示为：

式中：G(t)=diag[g(t−t1);g(t−t2);···;g(t−tn)]；x(t)=[x(t−t1),x(t−t2),···,x(t−tn)] ；Γ(t)为作用力指示向量，其作用是按轮距成反比的加权因子将轮对作用力分配到相邻的两个桥梁节点上；x(t)为平稳随机过程；G(t)为调制函数矩阵；n为轮对总数。

假定u˜i(t) 与u˜j(t)是系统任意时刻响应向量，则系统随机响应功率谱矩阵Suiuj(w,t)：

式中：hi(ti−τi,τi)为第i时刻脉冲响应函数。当i=j时，Suiuj(w,t)表示i时刻响应的自功率谱；当i≠j时，Suiuj(w,t)表示i、j时刻响应的互功率谱。

最后，根据随机响应功率谱密度矩阵，通过下式求得系统随机响应的均方根：

式中：Kz1、Cz1、Mw分别为车辆一系悬挂刚度、阻尼及轮对质量；Cc为轮轨间横向蠕滑而产生的附加阻尼；Sa、Sv、Sc分别代表轨道方向、高低、水平不平顺功率谱；b1为左右轨距中心距之半。横向力与竖向力标准差：

式中： σz、 σy分别代表竖向力与横向力的标准差；Szz、Syy分别表示由轨道高低不平顺引起的竖向力功率谱和轨道方向不平顺引起的横向力功率谱；Suz、Suy表示由轨道水平不平顺分别引起的竖向力功率谱和横向力功率谱；dw频率间隔带宽。

2.2 非平稳功率谱的谐波叠加法

采用谐波叠加法模拟随机过程是常用手段。土木工程领域常采用此方法模拟轨道、路面不平顺及随机风场。此方法要求功率谱函数是平稳过程[20]。对于轮轨力而言，由于桥梁位移变化具有时变特性，轮轨力功率谱也具有时变性，因此轮轨力功率谱函数，不能直接使用。

对于非平稳功率谱函数，引入调制函数，采用三角级数叠加法在每一时刻进行谐波叠加求其时域序列x(t)：

式中：x(t)为轮轨力时程序列；Sf(w)为轮轨力功率谱；g(t)调制函数；dw频率间隔带宽；wk为第k个对应频率值； ϕk为响应第k个频率的相位，一般按照0～ 2π间均匀分布取值； σF为非平稳激励下的标准差； σFs为平稳激励下的标准差。

2.3 极值理论

由概率论，服从极值分布需同时满足三个条件：观测对象是随机变量；这个随机变量的底分布保持不变，或者如果有任何变化，可通过数据变换减少其影响；观测到的极值是独立的[21]。对于列车脱轨系数分布而言，满足上述条件1 和条件3。对于高速铁路桥梁而言，桥梁挠跨比值受到严格控制，因此横、竖向轮轨力受轨道不平顺的影响十分显著，其分布高度近似服从高斯分布。脱轨系数是横向力与动轮重的比值，其安全指标为[−0.8，0.8]，对于动轮重和横向力，可以通过线性变换使其转化为标准正态分布，因此二者比值服从柯西分布。对于脱轨系数而言，满足极值分布的第二条件。由此可以得出，脱轨系数满足极值分布条件。

式中：u是位置参数； σ是尺度参数； γ为欧拉常数，取0.5772；Dx为极值分布的方差；Ex为极值分布的均值。

极值分布函数求解流程图，见图1。

图1 极值分布函数求解流程Fig. 1 Solution flow of extreme value distribution function

3 算例验证及计算结果分析

本文以秦沈客运专线现场实验为工程背景，通过与实验结果比较来验证本文模型与程序的准确性。因此，本文以秦沈客运专线狗河特大桥单跨24.6 m 标准简支梁为例；计算列车采用先锋号列车动车组，列车采用6 车编组，即动+拖+动+动+拖+动。运行时速270 km/h；线路横向偏心2.5 m，垂向偏心1 m，其轮轨力与桥梁作用位置关系图如图2 所示；轮轨间激扰采用由德国高速铁路低干扰谱变换出的轨道不平顺时域样本。车辆上桥前，先以与桥上相同的线路条件行驶150 m，待车体振动趋于稳定后进入桥跨结构，出桥后车辆再行驶450 m。桥梁与车辆参数见文献[19]。

图2 作用在桥梁上的力向量示意图 /cmFig. 2 Schematic diagram of force vector acting on bridge

由图3、图4 为本文计算结果和文献计算结果，二者比较可以看出，桥梁跨中竖向动位移幅值和变化趋势基本吻合。跨中竖向加速度变化趋势相近，幅值与文献计算结果有差异，主要原因是由于文献采用的是实测轨道不平顺，而本文采取德国高速铁路低干扰谱，轨道不平顺对桥梁加速度有较大影响，因此导致本文计算结果和文献结果稍有差异。通过上述计算结果与文献结果对比可验证模型与程序的正确性。

图3 跨中竖向动挠度Fig. 3 Vertical dynamic deflection in mid-span

图4 跨中竖向加速度Fig. 4 Vertical acceleration in mid-span

基于上述车-桥模型研究桥梁动力响应的随机特征。其中轨道不平顺为德国低干扰轨道不平顺谱，波长范围为1 m～80 m，相应的频率计算范围为0.9375 Hz～75 Hz，频域积分步长为0.9375 Hz，考虑轮对惯性力的影响。

3.1 桥梁变形对轮轨力平稳特性的影响

由轮轨密贴与横向线性蠕滑理论可知，竖向轮轨力含w2项，而横向轮轨力仅含w项，因此，高频对竖向轮轨力有较大影响。考虑轨道不平顺二阶导数项，结果如图5 所示。考虑轨道二阶导数项时，虚拟激励法求得的竖向力均方根值与100 次时域循环求得均方根值较为吻合。因此，采用虚拟激励法研究考虑轮对惯性力时的车桥系统随机振动，密贴模型能满足要求。

图5 首轮对竖向轮轨力均方根Fig. 5 Root mean square of vertical wheel rail force of first wheel pair

由图6 可知：横向轮轨力功率谱首次峰值出现在20 Hz 附近，频率越高，谱值越小；竖向轮轨力功率谱主要集中在中高频阶段，首次峰值出现在8 Hz 附近，当频率大于28 Hz，轮轨力谱值趋于稳定。对于密贴模型，轮轨力受高频影响较大。首轮对经过桥梁时间段为0.4 s～0.728 s，此段时间，横、竖向轮轨力功率谱波动不明显波动。

为探究桥梁位移变化对轮轨力功率谱平稳特性的影响，降低桥梁抗弯弹性模量，增大桥梁横、竖向变形。挠跨比取L/1500 与L/900，基于不同的挠跨比求动车首轮对轮轨力功率谱。由图6～图8 对比可知：横向力与竖向力功率谱幅值几乎无变化；增大桥梁变形对影响轮轨力功率谱值影响很小，桥梁位移变化对轮轨力功率谱平稳特性影响不显著，轮轨力功率谱平稳特性主要受轨道不平顺控制。

图6 动车首轮对轮轨力功率谱Fig. 6 Power spectrum of wheel rail force on first wheel of motor car

图8 桥梁挠跨比L/900 时的轮轨力功率谱Fig. 8 Wheel-rail force power spectrum of bridge with deflection span ratio of L/900

3.2 脱轨系数极值概率统计

车-桥系统动力响应安全检验通常采用3 倍标准差方法，此方法的前提条件是响应分布服从正态分布。脱轨系数服从柯西分布，因此采用3 倍标准差法则显然不合适。通过前节分析可知，脱轨系数也服从极值分布。本文通过轮轨力功率谱得到大量轮轨力时域样本，再结合轮轨力均值求得每次脱轨系数时域样本最大值，进而求得脱轨系数的极值分布。

采用谐波叠加法对轮轨力功率谱进行转换，得到10 000 组轮轨力时域样本。求得10 000 个脱轨系数最大值并进行极值分布检验。由图9 可以看出，动车、拖车首轮对脱轨系数最大值可以很好的服从极值I 型分布。置信度99.73%时，动车首轮对脱轨系数的极值为0.281；拖车首轮对脱轨系数的极值为0.385。

图9 首轮对脱轨系数极值分布Fig. 9 Extreme value distribution of derailment coefficient of first wheel pair

为了进一步验证极值I 型分布能很好的适用于脱轨系数最大值概率统计，抽取10 000 个样本中的一部分做极值I 型分布检验，从图10、图11 可以看出，极值I 型分布能够很好的适用与脱轨系数最大值概率统计。

图10 动车首轮对脱轨系数极值分布Fig. 10 Extreme value distribution of derailment coefficient of first wheel pair of motor car

图11 拖车首轮对脱轨系数极值分布Fig. 11 Extreme value distribution of derailment coefficient of first wheel pair of trail car

采用Monte Carlo 法验证，其计算样本数为1000。取每次计算结果的最大值，与极值I 型分布函数的99.73%置信度所对应极值对比。由图12可知，动车首轮对与拖车首轮对1000 次时程样本脱轨系数最大值各有2 次超过相应99.73%置信度对应极值。因此运用极值I 型分布求极大值是合理有效的。

图12 Monte Carlo 法与极值法所得脱轨系数最大值Fig. 12 Maximum value of derailment coefficient obtained by Monte Carlo method and extreme value method

3.3 基于不同车速的脱轨系数极值分布检验

列车运行速度对其脱轨系数影响显著。当波长范围一定时，速度越大，计算频率越高，横、竖向轮轨力越大。因此研究不同速度对脱轨系数的影响以及对极值I 型分布准确性的影响尤为重要。为研究车速变化对脱轨系数的影响，结合实际情况，选取270 km/h、300 km/h、330 km/h，图13 与图14 分别给出了极值I 型分布拟合与置信度99.73%的极值与100 次脱轨系数时程样本最大值关系图。

对不同车速下脱轨系数分布检验可以看出，列车速度越大，置信度99.73%的脱轨系数值越大，且对应的极值分布也越加发散，从图13(c)可以看出，当速度为330 km/h 时，置信度99.73%的脱轨系数极值为0.438。虽然发散程度存在差异性，但是脱轨系数分布均能很好的拟合于极值I 型分布。

图13 脱轨系数极值分布Fig. 13 Extreme value distribution of derailment coefficient

由图14 可看出，对于不同速度，100 次脱轨系数时程样本最大值均小于置信度99.73%的极值。因此极值I 型分布置信度99.73%所对应的极值可以较好的控制住不同速度下的脱轨系数最大值。

图14 Monte Carlo 法与极值法所得脱轨系数最大值Fig. 14 Maximum value of derailment coefficient obtained by Monte Carlo method and extreme value method

4 结论

本文采用虚拟激励法求得轮轨力非平稳功率谱，结合调制函数与谐波叠加法得到轮轨力时域样本，通过极值理论求得脱轨系数最大值的概率分布，并与Monte Carlo 法对比验证。算例结算表明本文方法的正确性，并得到以下结论；

(1)本文计算方法能合理地对不同速度下列车脱轨系数进行统计分析，采用极值I 型分布能够很好的适用于脱轨系数最大值的概率统计。

(2)桥梁变形对轮轨力功率谱平稳特性的影响很小，因此，对于高速铁路简支桥梁而言，轮轨力功率谱非平稳特性不显著。