基于DPM模型的高速列车转向架区域风雪流数值模拟

2023-09-25谢菲高广军何侃陈功

铁道科学与工程学报 2023年8期

谢菲，高广军，何侃，陈功

(1.湖南铁道职业技术学院，湖南株洲 412001；2.中南大学交通运输工程学院，湖南长沙 410075；3.中车株洲机车电力机车有限公司，湖南株洲 412001)

高速列车在寒冷的冬季运行时，轨道上的积雪会受到列车风影响，在列车的底部和尾部飞扬[1]。此外，列车转向架上的电机、制动钳夹等设备在运行时会产生热能。当空气中运动的雪粒子接触这些发热设备时，将会融化为水，进一步在高寒低温条件下凝结成冰。因此，列车在长期运行过程中，积雪会在转向架区域持续堆积，最终形成厚重的冰雪混合物，如图1所示。这一现象将增加转向架配重，危害转向架各关键部件的工作性能，进而对列车行车安全产生极大隐患。针对列车积雪结冰问题，JEMT[2]采用对环境无害的加热丙二醇作为材料加入除冰设备中，来改善瑞典等北欧国家的轨道积雪问题。PAULUKUHN等[3]针对俄罗斯冬季低温、多雪环境，为保证列车能够经得起无法避免的冰雪侵害，在除冰后的列车底部喷射水和丙烯类混合液，从而减少结冰。BETTEZ[4]针对日本铁路线的列车运行对高效除冰的要求，采用喷洒器作为除冰装置，有效解决了因积雪结冰而导致列车延误的问题。综上所述，来自于冬季处于高寒环境国家的研究人员提出了一系列的列车防积雪应对策略。然而，针对转向架区域的积雪结冰研究鲜有文献提及。高广军等[5-8]采用化学发热材料或物理结构改进等方式来改善转向架区域的积雪情况。丁叁叁等[9-11]均对列车转向架区域积雪结冰问题开展了流场数值模拟研究。然而，这些研究仅从流场进风量的角度间接预测转向架区域的积雪情况，并未考虑雪粒子在空气中的分布。雪粒子虽对流场运动具有一定跟随性，但气流出现扰动或者流速降低时，雪粒子的运动惯性将使得其对流场的跟随性减弱。在转向架区域，由于复杂几何模型影响，气流较为紊乱，存在大量湍流结构，且流速相对车身周围流场较低。与真实情况相比，仅从流场的角度考虑转向架区域积雪分布有失偏颇。因此，本文试图从风雪两相流的角度探索转向架区域积雪分布情况。关于风雪两相流，FUJII 等[12]针对日本铁路积雪情况进行研究，采用实验的方式分析了雪粒子运动的临界风速与雪密度及含水量之间的关系，提出雪粒子初始密度范围为36～132 kg/m3。TOMINAGA 等[13]采用数值仿真对雪粒子的堆积、侵蚀等情况进行模拟，并通过实验验证得出雪粒子的密度为100 kg/m3时数值仿真与实验结果具有良好的一致性。桑嘉宾[14]建立风雪流的数值模型，研究了混合粒径，得出雪粒子直径的影响因素，指出雪粒子直径一般分布在0～0.6 mm 之间，且主要集中在0.2 mm 附近。莫华美[15]针对屋面积雪分布建立了风雪流数值模型，利用该数值模型对屋面上的积雪分布规律进行数值模拟，指出雪粒子直径主要分布在0～0.4 mm 之间。TOMINAGA 等[16]采用CFD 数值模拟仿真的方法模拟雪粒子的堆积深度，并且采用风洞试验验证了模拟方法的可行性。BEYERS等[17]采用数值仿真模拟了雪粒子在建筑物周围逐步形成堆积的过程，有助于在建筑设计时对潜在的积雪情况进行处理。SERINE 等[18]将两维空间上的雪粒子数值模拟方法和风洞试验对比，并将其推广到三维空间上，从而确定雪粒子数值仿真的可行性。张强等[19]基于欧拉两相流理论和扩展传热模型，提出一种模拟飞机机翼积雪结冰的方法，并在实验中得到验证。GORDON 等[20-21]针对严寒天气的积雪现象，对雪粒子的形状、大小、速度及密度等进行研究，得出雪粒子的相关物理性质。SATO 等[22]针对在强风情况下的雪粒子破碎问题，采用风洞试验对积雪的密度进行研究，并在实地观测结果中验证其合理性。ZHOU 等[23]采用数值模拟技术对机场航站楼屋顶的雪粒子现象进行模拟，并用风洞试验对比验证数值模拟技术的可靠性，从而得出湍流模型、风速及风向等对积雪的影响。BEYERS等[24]采用数值方法模拟立方体周围的积雪情况，并通过实验进行对比验证。韩俊臣等[25]采用三维非定常雷诺时均模拟高速动车组转向架区域流场和雪粒子分布情况。娄振[26]采用准瞬态的方法模拟高速列车拖车转向架区域的积冰过程，研究转向架区域的积冰分布特征。上述文献表明，采用风雪两相流数值模拟仿方法可以再现积雪分布实际情况，具有较高的可行性。因此，本文采用数值模拟方法，复现列车转向架区域的流场流动情况，结合DPM 模型预测了雪粒子在空气中的运动轨迹，揭示了雪粒子在转向架上的分布情况，并从雪粒子随流场运动的角度描述了转向架区域积雪形成的原因，为高速列车转向架区域积雪结冰研究提供参考。

图1 转向架区域冰雪混合物堆积Fig.1 Accumulation of the ice and snow mixture in the bogie area

1 研究方法

1.1 物理模型、计算域及边界条件设置

本文基于转向架区域流场，建立了如图2所示的带有单个转向架区域的高速列车车体模型。考虑到实际运行中，头车的结构和形状与其他车不一致，不具有普遍性。因此，选择中间车型作为研究对象，在实践中能更好地估计一般转向架上的积雪结冰情况。该简化模型忽略了列车头、尾车等流线型结构对真实列车宏观流场产生影响的因素，进而着眼于转向架区域。车体的高度为3.6 m，宽度为3.2 m，长度为17.6 m。转向架区域距离车体迎风侧端面距离为6.9 m，使得车体底部流场在进入转向架区域之前得到充分发展，符合实际运行的列车转向架区域流场情况。同时，转向架区域距离车体背风侧端面同样为6.9 m，使得转向架区域的尾流得到充分发展，避免尾流对转向架区域流场产生干扰。此外，为真实模拟高速列车运行时车体底部的流场结构，增加了离地高度为0.2 m的钢轨模型。

图2 车体模型Fig.2 Geometry of the train

为模拟车体周围流场，建立了长60 m，宽20 m，高20 m 的计算区域，如图3 所示。车体上游处(平面ABCD)设置为速度入口，给定速度为55.56 m/s，使得雷诺数为1.325×107。车体下游处(平面EFGH)设置为压力出口，给定压力为0 Pa。车体两侧(平面ADHE和平面BCGF)及顶面(平面ABFE)均设置为固定壁面，模拟风洞环境。地面(平面DCGH)设置为滑移地面，给定与车速相同的滑移速度。车体及转向架壁面设置为固定壁面。车体距离气流入口15 m，车体两侧与计算域侧边界距离等宽。其中，模型的阻塞比约为2.8%，少于5%，满足一般钝体车辆计算需求[27]。计算时间步长设置为0.001 s，使得库朗数平均值小于1。

图3 计算区域Fig.3 Domain for simulation

1.2 网格设置

本文数值计算中采用六面体主导的网格。计算域采用开源计算流体力学分析软件OpenFoam 进行SnappyHex 网格离散，网格单元总数2 000 万。网格采用切割单元(cut-cell)方式，以级数递增方式从计算域边界向车体壁面逐层加密。车体壁面网格尺寸为0.02 m，转向架壁面网格尺寸为0.01 m。为准确模拟边界层流场，在车体和转向架壁面上布置了8 层附面层网格，近壁面第1 层网格高度为0.6 mm，使得壁面无量纲化尺寸y+处于30～300 之间。图4展示了转向架表面网格以及车身周围计算域剖面网格。

图4 数值计算网格Fig.4 Mesh for simulation

为了明确本文数值模拟不受网格密度的影响，基于相同的壁面法向y+，在流向和展向使用了不同的网格尺度，并采用了中等密度网格和精细密度网格进行对比。其中，中等网格和精细网格数量分别是2 000万和4 000万。如图5所示，为从车体底面中心线的2 个网格获得的时均压力系数Cp。在图中可以发现中等网格和精细网格的压力系数基本保持一致，说明本文采用的中等网格已经达到计算精度的收敛区间，具有较好的准确性。

图5 车体底面中心线压力系数对比Fig.5 Time-averaged Cp along the centerline of the bottom surface of the train body

1.3 流场数值模拟方法

针对高速列车周围流场，HE 等[28]指出采用大涡模拟(LES)能较好地捕捉湍流特性，但采用改进的延迟分离涡(IDDES)或者非稳态雷诺平均(URANS)方法同样能较为准确地模拟流场速度及压力分布。考虑到LES 或IDDES 对于网格分辨率要求较高，不适应于较高雷诺数条件下的数值仿真。此外，本文研究内容除了对流场运动方程求解之外，还需要额外求解雪粒子运动方程，计算量较大。因此，本文采用对网格尺度要求更低、在工程尺度得到广泛认可及应用的URANS 方法对流场部分进行求解。Realizablek-ɛ湍流模型能较好地捕捉流动分离，保持雷诺应力与真实湍流一致，使计算结果更符合真实性。因此，为适配壁面网格尺寸，选用了Realizablek-ɛ湍流模型，其控制方程见参考文献[29]。

1.4 雪粒子数值模拟方法

风吹雪是一种气固两相流，根据现代多相流计算理论，计算模型主要有欧拉两流体模型、混合模型以及拉格朗日粒子轨迹模型等[30]。拉格朗日粒子模型是另一种气固两相流计算方法，目前风吹雪数值模拟中应用较少。拉格朗日模型的优点是物理概念明确、结果直观，可计入粒子的轨道经历效应和历史效应，能得到粒子的详细运动信息[31]。因此，本文采用拉格朗日随机模型方法，对转向架区域的积雪分布进行数值模拟。

在本文中，通过对粒子的力平衡方程进行积分，来预测离散相粒子的运动轨迹。其中，作用在粒子上的外力主要包括气动力、重力、曳力等[32]。除此之外，其余力对粒子运动影响相对较小，在分析中可忽略不计[30]，其运动方程见式(1)：

对方程(1)积分就得到了颗粒轨道上每一个位置上的颗粒速度，见式(2)：

其中：x为粒子位移矢量，t为时间。粒子的轨迹是在离散的时间步长上逐步进行积分运算得到的，对式(2)进行积分，就可以得到粒子的轨迹。

对于车体、轮对等不易积雪部件，设定雪粒子边界为反射条件；对于构架、电机、端板、制动钳夹等易积雪部件，设定雪粒子边界为捕捉条件，认为转向架区域雪粒子全部被捕捉，以此来考虑积雪最严重情况。其中，雪粒子密度为100 kg/m3，雪粒子直径采用均匀分布为0.2 mm，雪粒子材料为water-liquid，列车运行速度为200 km/h。

此外，时间步长Δt由库朗数Co最大值控制，定义为Co=vΔt/Δx。其中，Δx是网格尺寸，v是网格内的流体速度。此外，隐式算法中的库朗数范围为0～100。如图6 所示，在本文数值模拟计算中，库朗数分布为0～11.85，且绝大部分集中在0和1 之间，基本满足计算要求。当t=2.00 s 时，气流多次通过该区域(风速55.56 m/s)，转向架中的风雪流场得到相对稳定和充分发展。并且，基于课题组多年计算转向架积雪的经验，2 s 的时间足以使得风夹带的雪粒子数量在转向架区域达到相对稳定。因此，本文的计算时间设置为2 s，时间步长为0.001 s，共计2 000步。

图6 库朗数等值线Fig.6 Contours of the courant number

2 数值模拟方法验证

关于本文中流场数值模拟方法的验证请参考文献[33]。本文仅针对数值模拟中雪粒子在转向架区域的分布情况进行验证分析，开展离散相风洞试验。该试验在中南大学高速列车研究中心的风洞实验室高速段进行。实验段的尺寸为3.4 m×1.0 m×0.8 m(长×宽×高)。由于真实雪粒子存在难制备、难保存等问题，试验中采用了与雪粒子密度相似的木屑粒子(密度约为100 kg/m3)作为替代物。木屑均经过筛网严格过滤，挑选出与雪粒子直径(0.2 mm)大小类似的粒子进行试验。试验段主要安装的装置包括高速列车原始转向架模型、高速列车原始车体模型、固定装置，其中固定装置又包括支撑杆、支撑座、分流板以及安装架等构成，如图7所示。

图7 模型的安装与固定Fig.7 Installation and fixation of the model

本次两相流风洞实验需要模拟雪粒子在转向架区域各关键部件的黏附情况，为此，在转向架和车体模型表面涂抹蜂蜜以黏附吹过的木屑。如图8所示，为涂刷蜂蜜后的转向架以及转向架各个细节部位的示意图。

图8 涂刷蜂蜜后的转向架Fig.8 Bogie painted with honey

为验证数值模拟方法的可靠性，且保障计算与试验的一致性，依据风洞试验室和试验模型的大小，建立了和风洞试验模型完全相同的数值模拟模型，且计算入口风速与风洞试验自由流风速一致。网格离散原则及数值模拟方法与本文第1节描述一致，由于简化了车体模型，网格总数为1 400万。

图9展示了风洞试验和仿真计算的转向架粒子堆积情况。从图9可以看出，风洞试验中木屑分布特点与数值仿真计算的积雪分布规律一致，粒子的堆积均产生在转向架各部件的迎风侧。因此，本文采用的数值仿真计算方法能较为准确地模拟雪粒子在转向架区域的分布情况。

图9 转向架底部表面积雪结果对比Fig.9 Comparison of snow results on the bottom surface of bogie

3 数值仿真结果分析

本节采用第1节所述流场及雪粒子数值模拟方法，针对转向架区域进行数值模拟仿真，得出转向架区域压力分布、流线分布、雪粒子分布及转向架表面积雪分布情况。

3.1 转向架区域压力分布

针对转向架区域主要部件的积雪结冰现象进行分析，在转向架区域作纵截面空间切片，如图10 所示，切片1 和4 在转向架轮对、制动钳夹位置，切片2和3在转向架电机位置。

图10 转向架区域切片位置Fig.10 Position of slices in the bogie area

如图11 所示，展示了转向架区域时均压力云图，来流方向为从左至右。由于切片1 和4 较为对称，为简略表达仅展示切片1，2 和3 上的时均压力分布。从图11 可以看出，轮对、制动钳夹、电机等突出部件在迎风侧由于受气流冲刷出现明显正压，而在这些部件的背风侧则为负压。转向架区域的压力分布体现了复杂几何结构对流场分布具有较大影响。因此，在进一步探讨雪粒子轨迹之前，需对流线分布进行分析。

图11 转向架区域压力分布Fig.11 Pressure distribution in the bogie area

3.2 转向架区域流线分布

图12 展示了转向架区域流线分布，来流方向为从左至右。从图12 可以看出，车体底部气流在转向架区域发生了偏转，部分气流在流经转向架前端时向上扬起进入转向架区域，并在转向架区域各个部件处形成多处漩涡，说明转向架区域流场具有较高湍流度。此外，气流进入转向架区域之后受复杂几何结构影响，速度明显降低，尤其在转向架轮对、制动钳夹以及电机等部位的流速接近为0。当流速降低时，雪粒子跟随性降低，当流速低于雪粒子起动风速时，重力占雪粒子所受合力比重增大。因此，当雪粒子随气流上扬进入转向架区域时，更易在低速区域附着至转向架各部件上形成堆积。值得注意的是，在转向架后端，部分气流在设备舱端板受阻，再次回流入转向架区域。

图12 转向架区域流线分布Fig.12 Streamline in the bogie area

3.3 转向架区域雪粒子分布

通过对转向架区域净风场流场情况的分析，已经初步判断出转向架区易形成雪粒子堆积区域。为进一步直观体现雪粒子在转向架区域的分布情况，图13 展示了雪粒子运动的快照，体现了雪粒子的运动历程。从图13 可以看出，在t=0.12 s 时刻，雪粒子进入到转向架区域的前部；在t=0.16 s时刻，雪粒子进入到转向架区域的中部，且集中在高度较低区域，对制动钳夹、电机等发热元件的冲刷不明显；在t=0.20 s 时刻，雪粒子经过整个转向架区域，更多的雪粒子进入到转向架的内部；当时间到达t=2.00 s 时刻，经过更长物理时间后，大量雪粒子进入转向架区域。此后，雪粒子分布趋于稳定。

图13 雪粒子运动时间序列快照Fig.13 Snap shots of the movement of snow particles

整体而言，在列车运行过程中，大量雪粒子跟随气流从列车车体底部流入，尽管大多数雪粒子到达并且通过了转向架区域，仍然有许多雪粒子随着上扬的空气进入了转向架区域。在转向架区域前方，雪粒子跟随气流运动，具有一定的跟随性；在转向架区域后方，雪粒子在运动惯性作用下和气流发生偏离。进入转向架区域的雪粒子流速降低、跟随性减弱，因此被转向架区域各部件捕捉。随着时间推移，雪粒子将会持续进入转向架区域，从而在各部件上产生堆积。

3.4 转向架表面积雪分布

如图14 所示，为展示转向架区域主要积雪部位，本文显示了转向架区域积雪浓度大于1×10-12kg/m3以上时的积雪覆盖情况，并采用速度渲染显示冲击到转向架区域时雪粒子的速度情况。图14 中，来流方向为从左至右。可以看出，在转向架底部区域产生大量积雪，且积雪更易堆积在转向架各部件迎风侧，而在各部件背风区域积雪相对较少。此外，积雪主要集中在电机、制动钳夹、轮对以及端板等部件底部迎风区域，而在转向架部件顶部和背风区域堆积相对较少。