冰雪条件下高速铁路道砟飞溅发生机理及影响研究

迟义浩，肖宏，张智海，郄录朝，张立军

(1. 北京交通大学 土木建筑工程学院，北京，100044；2. 北京交通大学 轨道工程北京市重点实验室，北京，100044；3. 中国铁道科学研究院集团有限公司铁道建筑研究所，北京，100081；4. 中国铁路兰州局集团有限公司运输部，甘肃 兰州，730000)

近年来，我国高速铁路实现了从无到有、从弱到强的跨越式发展，在综合交通运输体系中占有重要地位，极大地推动了我国经济的发展和社会的进步[1]。具有适应性强、成本低廉、易于维修等优势的有砟轨道[2]，在武九、长常、银西等设计速度为250 km/h 的高速铁路上均大量铺设，约占我国高速铁路总和的43%。在严寒地区，高速列车车体转向架位置常会出现冰块附着的现象，随外界温度的变化，在高速运行的动力荷载以及空气动力效应的耦合作用下，冰块脱落击打道床的情况容易发生，这会引起道床表面道砟发生飞溅现象，击打车体结构，对运营安全造成危害[3]。冰雪飞溅问题是我国高速有砟铁路在冰雪条件下遇到的关键问题，严重影响严寒地区高速铁路的运营安全。

近年来，国内外学者主要采用试验研究及数值模拟等手段，针对道砟飞溅的机理、影响因素及防治措施等进行了研究。在试验研究方面，IDO等[4-5]开展了道砟飞溅风洞及现场试验，主要研究了转向架附近位置的流场特性并开发了一种导流板装置，发现流场强度与车底光滑程度呈现负相关。KWON等[6]开展了道砟颗粒形状对道砟飞溅的影响研究，发现在相同条件下，道砟颗粒质量与表面积之比越大(如扁平状道砟颗粒)，对风载越敏感，发生飞溅的概率越高。SOPER等[7]基于现场试验，对轨道结构的流场特性以及振动特征进行了测试，他们发现不同列车类型对轨道动力响应特征影响较小，对流场特性影响较大。PREMOLI等[8]开展了有砟道床的实尺模型风洞试验，发现道砟飞溅数量与道床密实度呈负相关，与道床表面高度呈正相关，并指出轨道结构发生振动会加剧道砟飞溅现象。KAWASHIMA 等[9]采用空气炮进行冰雪飞溅试验研究，给出了飞溅的道砟数量、质量等与碰撞速度之间的关系。郄录朝等[10]在京沪高铁开展有砟轨道列车底部空气动力学现场试验，并结合仿真计算，对列车底部至轨道道床顶面之间的空气流动特性进行了分析。

在数值模拟方面，DENG 等[11]依托FRA 和AAR 提供的数据以及系统工程方法，利用可靠度进行冰雪飞溅机理量化分析，采用半定量分析方法来评估冰雪飞溅危害和相应限速标准。DING等[12]基于不连续变形方法(DDA 方法)对冰雪飞溅的机理进行了研究，发现随列车速度增大，激起的道砟颗粒数目和道砟颗粒运动的最大位移增大；入射角对道砟颗粒的运动方向的影响较大；冰块的形状影响道砟飞溅的数量和范围。LOPONEN等[13]采用基于黏附力的简化平衡方程，研究了列车底架位置冰块脱落所需的激励量，发现冰块脱落需要(20～2 000)g的加速度振幅。林建等[14]建立了列车风-轨枕-道床一体化耦合模型，根据仿真计算结果，提出了冰雪飞溅防治措施主要为线路结构优化、车体结构优化、站线除雪、接触网除冰。盖杰等[15]基于计算流体力学理论，构建了CRH5型高速列车模型，研究了冬季冰雪环境下运行工况和冰雪堆积对道床表面压力和道砟飞溅的影响，发现在交汇工况下，道床表面压力波动更为剧烈，且随冰雪层厚度增加，最大负压增大明显。

国内外学者对道砟飞溅开展了相关研究，但关于冰雪条件下道砟飞溅的研究较少，特别是对高速列车通过时，冰块脱落击打散体道床引起道砟飞溅力学行为和动态响应特征的建模理论及分析方法的研究更少，而揭示道砟飞溅运动轨迹及动态响应可以更好地制定防治措施。基于此，本文作者采用离散元法(DEM)、多体动力学(MBD)与计算流体力学(CFD)三者相耦合的方式，建立高速列车-有砟轨道冰雪飞溅精细化仿真模型，揭示高速列车运营条件下的有砟轨道流场分布特征，探究列车运行速度对冰雪飞溅行为的影响，探明列车风载及风振耦合下的冰雪飞溅动态响应特征，以期为严寒地区冰雪飞溅数值模型建立及防治提供参考。

1 高速列车-有砟轨道冰雪飞溅模型建立及验证

1.1 耦合仿真模型的建立

本文利用离散元法-多体动力学-计算流体力学(DEM-MBD-CFD)耦合分析的方法模拟道砟颗粒在冰雪作用下的运动及冰块脱落击打过程，精确考虑有砟道床的散粒体特性，以实现列车动力荷载的有效施加，真实还原高速列车-有砟轨道气动环境的流场特性。

1.1.1 有砟轨道离散元建模

根据TB 10621—2014《高速铁路设计规范》[16]确定有砟轨道模型尺寸，如表1所示。

表1 有砟轨道模型尺寸参数Table 1 Size parameters of ballast track model

有砟道床是有砟轨道结构的重要组成部分[17]，道砟颗粒的几何特征对道床的刚度、累积变形等力学特性具有显著影响，对道砟真实轮廓的还原重构是实现有砟道床精确模拟的关键[18-19]。为表征道砟粗糙表面及棱角特征，本文利用三维激光扫描还原了12 个具有典型特征的道砟颗粒的真实外形，并将球单元填充至重建的道砟颗粒三维廓形中，构建道砟颗粒的离散元模型。当道砟三维廓形的填充率达到90%以上时，再增加球单元的数量对仿真结果的影响程度很小[20]，故综合考虑计算效率与仿真精度，本文选用12～16个球单元黏结来模拟单个道砟颗粒，如图1所示。

图1 12个道砟颗粒模型Fig. 1 Models of 12 ballast particles

根据TB/T 2140－2018《铁路碎石道砟》[21]中的特级级配要求，绘制道砟颗粒级配曲线，如图2所示。表1所示为模型尺寸参数，采用“落雨法”建立有砟道床，利用SolidWorks 软件建立60 kg/m 钢轨和III 型混凝土轨枕几何仿真模型，并导入有砟道床模型中。综合考虑列车转向架轴距及计算效率，建立包含4根轨枕的有砟轨道模型，如图3所示。考虑边界效应的影响，主要对中间位置的2根轨枕及3个枕空区域进行分析。

图2 道砟颗粒级配曲线Fig. 2 Grading curves of ballast particles

图3 有砟轨道离散元模型Fig. 3 Discrete element model of ballast track

1.1.2 列车荷载施加

在高速列车荷载作用下，有砟轨道中不同位置轨枕对列车荷载有不同的分担比例。一般认为列车荷载由车轮下方的5根轨枕承担，其分担比依次为0.1:0.2:0.4:0.2:0.1[22-23]，可采用高斯函数对列车荷载进行拟合，同时考虑对多跨轨枕同时加载时存在荷载作用时间差现象，最终确定采用式(1)所示的多体动力学程序施加列车荷载函数p(v,t)：

式中：P0为列车轴重力；m和n分别为车厢数量和每节车厢的轮对数量；v为行车速度；t为时间；xij为车轮位置与车头的距离，考虑转向架轴距为2.5 m，车辆定距为17.375 m，车厢平均长度为25 m，选取列车车头位置为原点x0，故x11=2.562 5 m，x12=5.062 5 m，x13=19.937 5 m，x14=22.437 5 m；d为两跨轨枕间距，取为0.6 m；k为与模型第1根跨轨枕相隔的轨枕数。

1.1.3 流场分析模型建立

本文以银西高铁上主要运行的CRH380B 型动车组列车为例，根据CRH380B 型列车三视图及相关参数，采用SolidWorks 三维建模软件建立列车模型。由于转向架位置易存在冰雪附着现象，因此，需重点关注轨道与车底之间的流场特征，综合考虑仿真效果及计算效率，采用头车-中间车-尾车的3节编组形式[24-25]，精确模拟列车底部结构如转向架整体结构，忽略列车表面存在的受电弓、车窗等多处凸起[26]及转向架的弹簧、垫片等细部结构，建立高速列车几何模型，如图4 所示。图4中，列车总长度为76.525 m。

图4 高速列车几何模型Fig. 4 Geometric model of high-speed train

采用基于滑移网格技术的分析方法模拟列车运行，建立高速列车移动计算域和有砟轨道固定计算域，实现对轨道流场的真实还原。滑移网格的基本原理为通过在2 个区域的交界面位置设置interface 单元，实现滑移运动中2 个区域之间的数据传递，通过滑移网格的运动，实现列车与轨道结构之间的相对运动。

考虑到计算量、流场充分发展和模拟的真实性，经过试算得到计算域尺寸，其中高速列车移动计算域长度为400 m，宽度为5 m，高度为5 m。有砟轨道固定计算域长度为200 m，宽度为60 m，高度为30 m。初始时刻，列车头车鼻尖距离轨道水平方向距离为20 m，沿轨道线路中心线正方向行驶。高速列车-有砟轨道流场分析计算模型如图5所示。

图5 高速列车-有砟轨道流场分析模型Fig. 5 Field analysis model of high-speed train-ballast track flow

利用ICEM软件进行模型的网格划分时，为保证计算精度同时提高计算效率，在靠近列车及轨道结构位置采用八叉树算法划分四面体非结构化网格，在远离列车的流场区域采用六面体结构化网格，通过共节点设置进行不同网格间的数据交换，同时对车体表面、钢轨、轨枕等轨道结构进行6层边界层设置，生成的网格模型如图6所示。

图6 网格模型Fig. 6 Model of grid

考虑到多面体网格可以在保证精度的情况下显著降低网格数量，大幅提高计算效率[27-28]，因而，利用ICEM导出生成的网格，在Fluent里转换为多面体网格进行求解计算，最终共生成约3 000万个网格单元。

将所有入口和出口边界(即高速列车移动计算域和有砟轨道固定计算域的前后截面)均设置为压力出口边界；有砟轨道固定计算域的流场顶面，左右侧界面均设置为对称边界；将各计算子域之间的滑移网格交界面设置为interface边界；混合网格交界面设置为内部面interior。此外，高速列车模型表面以及钢轨、轨枕、道床等轨道结构表面均设置为无滑移壁面边界。

本文列车最大运行速度为350 km/h，马赫数小于0.3，故认为空气是不可压缩的，采用RNGk-ε湍流模型，利用SIMPLEC 压力速度耦合算法以提高计算效率和精度[29-30]，设定计算时间步长为0.005 s，时步内迭代次数为30次。

1.1.4 冰块击打过程模拟

采用CFD-DEM单向耦合以及DEM-MBD双向耦合的方式进行冰块击打过程的模拟计算，具体模拟过程如下。

1) 采用计算流体力学软件，进行高速列车-有砟轨道流场分析模型的瞬态计算，在计算完成后，编写UDF 程序，提取流场中各个网格节点位置的速度。由于流场分析模型的长度远比有砟轨道离散元模型的长度大，且在每个固定时步内均会产生流场数据文件，即头车前转向架附近区域的轨道流场强度及波动特征较为显著[31-32]，利用Matlab编制函数，只提取该区域内的轨道流场数据。

2) 采用API 编译“场力模型”，导入过程1)中提取的轨道流场数据，对有砟道床道砟颗粒单元及冰雪颗粒施加“场力”作用，与此同时，基于DEM-MBD 耦合分析方法，将列车荷载直接施加在轨枕上，依托“wall”单元实现信息传递。

3) 基于现场拍摄的车厢位置处的冰块形状，建立冰块离散元模型，冰块质量约为0.3 kg，如图7所示，将其导入有砟轨道离散元模型中。假设冰块从中间车第一转向架后的车厢中部位置上掉落，在预设位置建立正方形颗粒工厂，采用cubic方法设置冰块的精确位置，并给定冰块初速度与列车运行速度相同，最终建立基于DEM-MBDCFD耦合的高速列车-有砟轨道冰雪飞溅模型，如图8 所示，图中蓝色区域表示轨道流场速度分布范围。

图7 冰块离散元模型的建立Fig.7 Establishment of models of ice block discrete element

图8 基于DEM-MBD-CFD耦合的冰雪飞溅模型Fig. 8 Ice and snow splash model based on DEM-MBD-CFD coupling

1.2 接触模型及参数

颗粒与颗粒之间、颗粒与几何体之间的接触均采用Hertz-Mindlin 接触模型，可由式(2)和式(3)计算两者之间的相互作用力，包括法向接触力Fn和切向接触力Ft：

式中：Ea、νa、Ga和Ra分别为颗粒与颗粒间接触的弹性模量、泊松比、剪切模量以及半径；Eb、νb、Gb和Rb分别为颗粒与几何体间接触的弹性模量、泊松比、剪切模量以及半径；α为法向重叠量；δ为切向重叠量。

冰块在风场影响下受到的曳力F的计算公式为

式中：Cdrag为曳力系数；Ap为颗粒投影面积；ρf为流体密度；vrel为相对速度。曳力系数计算公式与雷诺数Re有关：

式中：d为特征长度；μ为流体动力黏度。

曳力系数的计算公式如下。

1) 当Re≤0.5时，

2) 当0.5

3) 当Re>1 000时，Cdrag=0.44。

模型参数主要包括基本物理属性参数以及接触参数，材料基本属性主要通过参考文献[2，17，19-20，33]进行拟定；为保证接触参数的可靠性，参考文献[12，34]中数据并通过控制变量法进行大量试算后拟定，耦合仿真模型在计算过程中所涉及的主要参数如表2所示。

表2 模型参数Table 2 Model parameters

1.3 模型验证

为验证耦合模型的准确性，采用风压传感器在银西高铁实测列车通过时道床表面的空气动力效应，现场测试照片如图9所示。在提取冰雪条件下，当列车运行速度不同时道床表面中心测点(道心)正压及负压的实测数据，绘制试验与仿真条件下的道心正压、负压随车速的变化规律，如图10所示。

图9 现场测试照片Fig. 9 Photoes of testing on spot

图10 风压现场试验值与仿真值对比Fig. 10 Comparisons of field test and simulation results of wind pressure

由图10 可知：道心正压、负压随车速变化的仿真值与试验值差距不大，道心正压试验值和仿真值相对误差平均为5%左右，道床负压试验值和仿真值相对误差平均为10%左右，误差较小。从变化规律可知，道床正压、负压试验值与仿真值随车速变化规律大致相同。

为进一步对模型的接触参数进行标定，验证冰块击打过程的可靠性，参考KAWASHIMA 等[9]的空气炮现场试验，建立颗粒床离散元模型，模拟一个同等大小的速度为100 km/h 的长×宽为600 mm×85 mm的冰块的碰撞行为，并将碰撞击打室内试验计算结果与试验观测结果进行比较，如图11所示。

图11 碰撞击打室内试验值与仿真值对比Fig. 11 Comparisons of indoor test and simulation results of collision striking

由图11 可知：碰撞击打室内试验离散元仿真值较好地显示出与试验观察到的相似的道砟飞溅现象，验证了道砟颗粒与冰块碰撞时的动态力学行为。综上可知，在误差允许的范围内，本文建立的模型是可靠的，可用于后续的冰块击打引起道砟飞溅的计算分析。

2 计算结果分析

2.1 风压特征分析

当高速列车通过有砟轨道模型正上方时，高速铁路-有砟轨道流场风压云图如图12所示。

图12 车体及轨道结构表面风压分布云图Fig. 12 Cloud maps of wind pressure distribution on surface of train body and track structure

由图12 可知：风压主要呈现出“正压—负压—负压—正压”的分布特征，车头前方存在以车头前端为圆心的较大的圆形正压区域。车头后方存在明显的负压区域，大致位于前转向架附近位置。车尾与车头大致呈对称分布，但其幅值明显比车头的低，且分布范围较小。从整体看，在列车中部转向架区域和车厢间隙连接处也存在一定的风压变化，这些结构不平顺部位容易引起流场急剧变化，形成复杂的湍流运动。在严寒地区，风中携带的雪花颗粒易受气流扰动的影响，流入并逐渐附着在转向架结构表面，进而当温度发生改变时，在列车动力荷载和空气动力效应的耦合作用下，列车高速通过时易引起冰雪飞溅问题。

为探究轨道流场特征随列车通过的变化情况，以列车运行速度为350 km/h 为例，列车运行长度总长为296.525 m，运行时间约为3.05 s，提取道床表面中心测点的风压，得到风压时程曲线，如图13所示。

图13 道床表面中心测点风压时程曲线Fig. 13 Time-history curves of wind pressure at the central measuring point on surface of ballast bed

由图13 可知：在列车运行过程中，道床表面中心存在2 个较为显著的压力波动：

第一个主要压力波动产生在车头前端区域，最大正压为2 004.5 Pa，而后方存在正负压交替的现象，压力梯度非常大，最大负压产生于头车前转向架位置附近，为-1 520.4 Pa。这是由于空气在流动时，受到车体的阻滞作用而产生的高速气流冲向道床和车体与道床间隙中引起的。

第二个主要压力波动产生于车尾后端区域，当车尾后转向架通过监测位置时，道床表面出现显著的负压波动，当车尾鼻端通过监测位置时，道床表面出现显著的正压波动，最大正压为1 096.2 Pa，最大负压为-1 059.9 Pa。因此，在对高速列车-有砟轨道冰雪飞溅模型进行计算分析时，应重点关注头车前转向架位置的流场速度。

2.2 不同行车速度影响研究

本节仅考虑在列车空气动力荷载作用下，冰块脱落发生击打道床造成道砟飞溅的行为，探究不同行车速度下列车风载对冰雪飞溅行为的影响。列车运行速度从150 km/h 递增至350 km/h，以50 km/h递增。

根据现场测量结果，CRH380B 型车的车厢底板距离钢轨顶面的高度为275 mm，设置道砟飞溅的监测区域为长方体，区域长×宽×高为2 600 mm×1 500 mm×275 mm，包含中间2 根轨枕和3 个枕空区域，监测区域底面与钢轨顶面平齐，如图14所示。

图14 道砟飞溅监测区域Fig. 14 Monitoring area of ballast flight

提取不同列车运行速度下，冰块脱落击打道床激起道砟飞溅的数量，其随时间的变化曲线如图15所示。

图15 不同速度下道砟飞溅时程曲线Fig. 15 Time-history curves of ballast flight at different speeds

由图15 可知：当掉落的冰块进入监测区域与道砟发生碰撞击打时，冰块会迅速发生反弹，并使与之接触的道砟颗粒发生飞溅，道砟颗粒在与冰块碰撞后获得初速度，在碰撞作用及列车风载影响下运动。

图16 所示为列车风载下，冰雪飞溅道砟数量与行车速度的关系。由图16 可知：列车运行速度是诱发冰雪条件下道砟飞溅的主要影响因素之一；随列车运行速度提高，飞溅道砟数量显著增加；道砟飞溅击打列车底板的临界速度为250 km/h，当行车速度从150 km/h 提升至350 km/h 时，飞溅道砟数量从1个最多增至16个；当行车速度从250 km/h提升至350 km/h时，击打列车底板的道砟数量从3个最多增至6个。

图16 列车风载下冰雪飞溅道砟数量与行车速度的关系Fig. 16 Relationship between number of ballast in snow flight and driving speed with train wind load

进一步提取不同行车速度下，列车底板受到的冲击力，如表3所示。

表3 列车底板受到的冲击力Table 3 Impact force on train floor

由表3可知：当行车速度为150 km/h和200 km/h时，未发生因道砟飞溅造成列车底板受击打现象；当行车速度为250 km/h 时，冲击力为2.80 kN，而当行车速度为300 km/h时，冲击力增大为7.22 kN，增大了4.42 kN，增长约157.9%；当行车速度为350 km/h时，冲击力增大为14.60 kN，较300 km/h时增大了7.38 kN，增长约102.2%，由此可看出，冲击力受行车速度的影响极大，当行车速度超过250 km/h时，冲击力随车速增长而成倍增加。

2.3 风振耦合影响研究

本节综合考虑列车风载和列车动力荷载的共同作用，结合CRH380B 型车运行实际，对轨枕施加轴质量为17 t，编组为8 的列车荷载函数。绘制不同列车运行速度时风振耦合条件下的冰雪飞溅道砟数量与行车速度的关系，如图17所示。

图17 风振耦合下冰雪飞溅道砟数量与行车速度的关系Fig. 17 Relationship between number of ballast in snow flight and driving speed with wind-vibration coupling

由图17 可知：在风振耦合作用下，飞溅道砟数量、道砟飞溅最大位移以及击打列车底板的道砟数量均明显提升；当行车速度达到200 km/h时，冰块击打引起的道砟就能够击打到列车底板，究其原因是在风振耦合作用下，道床振动加剧，使得道床表层道砟颗粒之间接触点的数量减少，导致表层道砟颗粒之间排布变松，稳定性减弱，冰块击打下飞溅的道砟数量增加，且一旦道砟颗粒飞离道床，在列车风载作用下，颗粒会在负压的作用下进一步加速吸起，击打车厢底板结构，影响行车安全。

对比列车风载和风振耦合条件时，在不同行车速度下，列车底板受到的冲击力关系如图18所示。

图18 列车底板受到的冲击力对比Fig. 18 Comparison of impact force on train floor

由图18 可知：当仅考虑列车风载作用时，冰雪飞溅击打列车底板的临界速度为250 km/h，随行车速度不断增大，击打时对列车底板的冲击力成倍数增大。在风振耦合共同作用下，临界速度降低为200 km/h，且冲击力随行车速度增大呈二次多项式分布递增，可由式(9)拟合得到。最大冲击力约为仅考虑列车风载情况下的2.2倍。从两者间的对比可看出，列车动力荷载可加剧冰雪飞溅程度。

式中：y为列车底板受到的冲击力；x为列车运行速度。

3 结论

1) 以CRH380B 型车为例，通过DEM-MBDCFD耦合分析方法，构建高速列车-有砟轨道冰雪飞溅精细化仿真模型，实现了冰块脱落击打道床，引起道砟飞溅过程的有效模拟，并通过现场实测数据，验证了模型的准确性。

2) 在冰雪条件下高速列车运行过程中，道床表面中心风压呈“正压—负压—负压—正压”的波动特征，且头车引起的正压和负压幅值均显著比尾车的幅值高；转向架位置处负压较大，流场波动明显，易使得冰雪在气流扰动的影响下，积聚附着在转向架结构表面。

3) 列车运行速度是诱发冰雪条件下道砟飞溅的主要影响因素之一。飞溅道砟数量与列车运行速度呈正相关，在列车风载作用下，当行车速度从150 km/h提升至350 km/h时，道砟飞溅数量从1个增长至最多16个，且基本呈线性增加关系。

4) 在风振耦合作用下，道砟飞溅击打到车厢底板的临界行车速度从仅在列车风载作用下的250 km/h降低到200 km/h，飞溅道砟总数量增长了约56.5%，列车底板受到的最大冲击力提升了118.6%，列车荷载显著加剧冰雪飞溅现象。因此，建议，在冰雪严重的高速铁路线上，为减少道砟飞溅击打车体底部，在不采取任何防治冰雪飞溅措施的条件下，运营速度应低于200 km/h。