基于DEM-MFBD-CFD方法的冰雪飞溅击打应答器力学研究

迟义浩, 肖 宏, 郄录朝, 张智海, 路 锟

(1. 北京交通大学 土木建筑工程学院,北京 100044; 2. 北京交通大学 轨道工程北京市重点实验室,北京 100044;3. 中国铁道科学研究院集团有限公司铁道建筑研究所,北京 100081;4. 中国铁路兰州局集团有限公司兰州车辆段,兰州 730050)

高速铁路以其高效性、便捷性、安全性的突出优势,在综合交通运输体系中占有重要地位[1]。我国幅员辽阔,地区差异性大,随着“交通强国”战略的日益推进,高速铁路在我国东北、新疆北部等寒冷地区也逐渐建设发展。但高寒地区由于冬季独特的低温冰雪气候,在列车开通运营后,会出现积雪、融冰随快速运营的列车掉落击打道床及应答器等,严重影响高速铁路运行效率及行车安全[2-3]。

国内外学者针对冰雪飞溅问题的产生机理、防治措施等已运用现场试验和数值模拟等手段开展了研究。Saussine等[4-5]开展了轨道结构的风洞试验,通过观察钢轨内侧枕间区域处的道砟颗粒运动情况,分析研究风场对道砟颗粒的运动影响;此外,基于随机过程法进行了冰雪飞溅的可靠度及相关防护研究,进而提出了冰雪线路运行限速及轨道结构设计方案等[6-7]。Ido等[8]将人工雪颗粒引入到有砟轨道风洞试验中,对线路积雪及自然降雪的现象进行了有效模拟,并开发了一种导流板装置,结果表明,车厢底部转向架区域的积雪情况在采用该装置后得到了显著的缓解。北欧国家[9-10]进行冰雪线路的系统研究,包括接触网除冰、冰雪堆积与除雪、道岔融雪,编制了冬季铁路运营与养护维修规范,重点对防雪、除雪、融冰和道床断面选型进行综合防治。林建等[11-13]通过对高速列车和轨道线路进行结构优化,如限制道砟颗粒粒径,取消砟肩堆高,应用聚氨酯固化道床等措施来控制严寒地区因冰雪问题引起的道砟飞溅现象。冯永华等[14]采用流体动力学仿真,分析了转向架区域的速度场与压力场,结果表明:随着速度增加,转向架区域中夹带的雪粒子增加,相应部位的积雪、结冰问题更为严重。余卫巍等[15]分析冰雪天气下动车组车底部掉落冰块,击中在线路中心安装的应答器和电源线,将造成应答器被击碎和电源线变形、断线等现象。

综上可知,针对冰雪飞溅问题,国内外学者已开展了研究,但主要关注的是冰雪飞溅对轨道结构和车体等的影响,目前关于冰雪飞溅问题对应答器的影响方面的文献报道较少,缺乏对列车高速运行下冰块脱落击打应答器的量化认识,无法科学指导养护维修。根据目前我国已经开通的哈大高铁、银西高铁等高速铁路线路的现场实际,已出现多处冰块脱落击打应答器的现象,甚至严重损坏了这些设备,部分地段不得不采取降速运营的措施,大大影响了列车的运营效率与安全性。基于此,基于离散元-多柔性体动力学-计算流体力学(discrete element method-multi flexible body dynamics-computational fluid dynamics,DEM-MFBD-CFD)三者耦合的分析方法,建立车厢底板结冰脱落击打应答器的仿真模型,探究行车速度、横风条件、冰块质量等因素对应答器击打的影响,提出防止结冰击打应答器的预防性措施,以期为严寒地区冰雪飞溅影响应答器的养护维修提供一定的依据。

1 DEM-MFBD-CFD耦合分析方法

冰雪飞溅击打应答器过程,需要考虑高速列车与有砟轨道间的相互作用,涉及流体力学、颗粒物质理论等多学科,过程极其复杂。为充分反映、更好模拟这一复杂过程,采用DEM-MFBD-CFD联合仿真的分析法,通过有砟轨道DEM模块作为中间媒介,依托DEM-CFD耦合,实现流场与道床之间的相互作用,通过DEM-MFBD耦合,实现列车荷载下有砟道床与轨枕的相互作用,并真实反映应答器的受力情况。依托耦合并行算法使得3个模块之间的动量、质量、受力等信息进行数据传递,三者间的耦合流程如图1所示。

(1)采用DEM建立有砟道床结构

将研究介质有砟道床划分为一系列离散单元,即不同的道砟颗粒,并赋予每一个道砟颗粒单元独立的几何物理特性,利用颗粒单元间的接触本构模型和接触关系,依托牛顿第二定律在每一个时间步长内计算运动情况,进而更新所有颗粒的位置信息,完成迭代运算。

(2)采用MFBD建立轨道应答器结构

应用有限元网格划分技术,对应答器结构进行节点划分及单元化处理,仿真过程中通过不断迭代计算结构矩阵,进而实现针对结构应力应变信息的精准分析。轨枕由于自身变形相对于结构运动位移可以忽略不计,采用刚体模型。

(3)采用CFD模拟高速列车运行产生的风场

基于有限体积法,将质量、动量、能量三大基本守恒定律的积分型方程或偏微分型方程转换成离散代数形式,然后求解这些方程得到离散的空间和时间点上的流场数值,采用滑移网格分析法,以准确模拟高速列车与有砟轨道之间的相对运动。

2 模型建立及验证

2.1 车厢底板结冰脱落击打应答器模型建立

2.1.1 有砟轨道道床离散元模型

既有研究表明,道砟颗粒的棱角和形态特征是准确建模的关键参数[16-18],为更好地反映道砟颗粒间咬合作用,保障有砟道床结构的真实性、合理性。基于三维激光扫描技术,准确获取了道砟颗粒的不规则外形轮廓,通过球状颗粒堆积法生成不规则道砟颗粒的模板,如图2所示。由图2可知,随着球颗粒数的增加,道砟颗粒模板更接近真实的道砟颗粒形状,但同时也存在着计算耗时高、占用内存大等局限性,综合考虑道砟颗粒的模拟精度和计算效率,结合国内外研究经验[19-21],在实际建模过程中利用9个球体填充形成道砟颗粒簇模拟一个道砟颗粒。

图2 单个道砟颗粒建模Fig.2 Single ballast particle modeling

根据TB/T 2140-2008《铁路碎石道砟》[22]规范要求,建立特级道砟颗粒级配,颗粒级配曲线如图3所示。选取10种具有典型颗粒特征的道砟样本。

图3 颗粒级配曲线Fig.3 Particle grading curve

根据TB 10621-2014《高速铁路设计规范》[23]建立高速铁路有砟轨道结构离散元模型,轨枕为III型混凝土轨枕,60 kg/m钢轨,道床顶面低于承轨面40 mm。为准确模拟有砟道床,根据图3的道砟颗粒分级级配,基于落雨法生成道砟颗粒,综合考虑计算效率和计算分析需求,采用分层压实法建立包含4根轨枕的有砟道床模型,如图4所示。道床的顶面宽度为3 600 mm,厚度为350 mm,边坡坡度为1∶1.75,轨枕长度为2 600 mm,间距为600 mm,模型中共包含107 692个道砟颗粒簇。

2.1.2 应答器柔性体模型

建模时,首先参照应答器尺寸图,采用三维绘图软件建立应答器三维模型图,接着将三维几何模型以“.step”格式导入到MFBD模块中,对应答器材料的基本物理属性进行设置,采用8节点六面体Solid8单元进行离散,最小网格单元尺寸约为0.1 mm,划分得到的网格数目为10 864个,建立柔性应答器模型,如图5所示。

图5 应答器柔性体建立(mm)Fig.5 Establishment of transponder flexible body(mm)

将应答器按照现场实际情况固定在轨枕上,以“.wall”文件的形式导入到DEM模块中,实现计算过程中数据的传输,最终建立起DEM-MFBD耦合分析模型如图6所示。其中模型中主要部件的基本物理力学参数如表1所示[24-25]。

表1 基本物理力学参数Tab.1 Basic physical mechanics parameters

图6 DEM-MFBD耦合分析模型Fig.6 DEM-MFBD coupling analysis model

2.1.3 流场分析模型

(1)控制方程

在本文涉及的流场分析中,设置的风速最大值为350 km/h,高速列车外流场运动马赫数均小于0.3,故可以认为流体是非压缩的,选用不可压缩流体的k-ε两方程湍流模型数值分析方法分析高速铁路-有砟轨道表面流场空气动力学特性,遵循密度为常数的定常流动连续性方程如式(1)所示[26]

∮ρ(V·n)dA=0

(1)

式中:V为流体速度矢量,m/s;n为法线方向;A为控制体面积,m2。

动量方程如式(2)所示:

∮ρV(V·n)dA=∮ρfdτ+∮pndA

(2)

式中:f为体积力,N/m2;τ为微元体;pn为作用于微元体表面的压强,Pa。

(2)几何模型

我国高速铁路里程长,车型种类较多,不同的车型其外观线型也存在差异,以高寒地区常运行的CRH380B型车为例建模。考虑到高速列车中间车体的截面形状相同、气动力的变化基本一致,为提高计算效率,取一节中间车采用“头车-中间车-尾车”3节编组(25.850 m+24.825 m+25.850 m)的形式[27],建立的仿真模型如图7所示。

图7 高速列车仿真模型Fig.7 Simulation model of high speed train

采用基于滑移网格技术的分析法,建立高速列车移动计算域和轨道结构固定计算域,通过设置交界面进行不同计算域之间的数据传输。考虑到计算量、流场充分发展和模拟的真实性,经过大量试算,确定计算域尺寸。高速列车移动计算域长度为400 m,宽度为5 m,高度为5 m,轨道结构固定计算域长度为200 m,宽度为80 m,高度为50 m。本模型采用混合网格方式划分网格,共计生成约3 000 万个网格单元,网格划分结果如图8所示。

图8 流场分析模型及网格划分(m)Fig.8 Flow field analysis model and meshing(m)

采用RNGk-ε湍流模型计算流场参数的变化情况,无横风时,设置前后截面为压力出口边界,轨道结构固定计算域的顶面及左右侧界面为对称边界,交界面设置为interface边界条件建立边界关联,其他结构表面设置为无滑移壁面边界条件。有横风时,将轨道结构固定计算域的左右侧界面更改为速度入口和压力出口边界。

利用SIMPLEC压力速度耦合算法以提高计算效率和精度。设定计算时间步长为0.005 s,设置时步内迭代次数为30次,使在这一时步内的流场可达到收敛状态。

2.1.4 击打过程模拟

基于DEM-MFBD-CFD耦合算法,结合现场实际,假定冰块从中间车第一转向架后的车厢中部位置上脱落,同时进行高速列车-有砟轨道流场分析计算与施加列车动荷载作用,利用MATLAB软件编制函数,提取列车经过时的流场数据,实现冰块在列车运行及风场作用下脱落击打应答器这一过程的有效模拟,从而精确反映在空气动力荷载、轨道动力荷载作用下的击打运动行为,最终建立的模型如图9所示。

图9 基于DEM-MFBD-CFD的车厢底板结冰脱落击打应答器模型Fig.9 DEM-MFBD-CFD based model of carriage floor ice-shedding attacking trackside equipment

2.2 模型验证

接触模型采用Hertz-Mindlin(no slip)模型,参照文献[28-29],模型中各部件的参数取值,如表2所示。

表2 模型参数Tab.2 Model parameters

参照风洞试验结果[30],开展行车速度为350 km/h时的高速列车-有砟轨道流场仿真计算,提取道床上方不同垂向位置和横向位置处的风速值,绘制仿真与试验结果对比曲线,如图10所示。

图10 风速结果验证Fig.10 Verification of wind speed results

由图10可知,从线型变化规律来看,模型计算得到的道床上方垂向风速分布特征和横向风速分布特征均与风洞试验结果得到的规律较为一致。从数值上看,差值也较小,最大误差约为4.0%。

为进一步验证模型的正确性,在现场布设风压传感器,测量不同列车运行速度下的道床表面中心测点正压数据,绘制仿真与试验对比结果如图11所示。

图11 风压结果验证Fig.11 Verification of wind pressure results

由图11可知,不同列车运行速度下,试验值与仿真值的差值在5.0%以内,两者较为接近,差值较小,说明模型在不同车速下的风压结果是可靠的。

综合图10、图11可知,从风压和风速两个角度都验证了所建立的模型的正确性,可以用于后续分析计算。

3 计算结果分析

3.1 行车速度影响分析

综合考虑我国列车运行速度实际情况,本节针对150 km/h,200 km/h,250 km/h,300 km/h,350 km/h不同时速,计算分析在列车空气动力效应以及列车动荷载的风振耦合作用下,车厢底板结冰脱落击打应答器的情况。

以350 km/h的列车行车速度,3 kg的冰块质量为例,绘制车厢底板结冰脱落击打应答器前后的应答器表面应力分布云图,如图12所示。

图12 击打前后应答器表面应力分布Fig.12 Surface stress distribution of transponder before and after striking

绘制击打前后应答器表面击打位置处的应力变化时程曲线如图13所示。

图13 应力变化时程曲线Fig.13 Time history curve of stress change

由图12、图13可知,冰块击打为点接触,应力呈现从中心向四周扩散的趋势,最大应力值约16 MPa。

绘制从150 km/h递增到350 km/h的不同行车速度下,应答器表面击打位置的应力分布时程曲线如图14所示。应答器表面的最大应力随行车速度的变化曲线如图15所示。

图14 不同行车速度下应答器表面击打位置应力分布时程曲线Fig.14 Time-history curve of stress distribution on beat position of transponder surface under different driving speeds

图15 最大应力与行车速度的关系Fig.15 Relationship between maximum stress and driving speed

由图14可知,不同行车速度下,冰块对应答器的击打应力显著不同,速度越快,击打应力越大。在行车速度为150 km/h的情况下,应答器表面所受到的最大应力数值为4.487 MPa,当行车速度增大到350 km/h时,最大应力增大为15.591 MPa,增长了11.104 MPa,约为150 km/h时的3.5 倍。这表明行车速度对应答器的击打具有主要影响,建议在冰雪严重地段,可采取降速措施,以更好地降低冰雪击打的影响。采用幂函数对不同行车速度下,应答器受击打的最大应力结果进行拟合,关系式如式(3)所示,决定系数R2值为0.991,拟合效果良好(见图15)。

y=0.002 4x1.51

(3)

式中:y为最大应力,MPa;x为列车运行速度,取值范围为150～350,km/h。

3.2 横风条件影响研究

考虑到我国目前高速有砟铁路的最高运营速度是250 km/h,下面分析行车速度为250 km/h且考虑在5 m/s,8 m/s,11 m/s,14 m/s,17 m/s,20 m/s的横风风场影响、列车空气动力效应以及动荷载的风振耦合作用下,车厢底板3 kg的结冰脱落击打应答器的情况,计算其表面击打位置处的应力大小变化的时程曲线,如图16所示。

图16 不同横风风速下应答器表面击打位置应力变化时程曲线Fig.16 Time-history curve of stress change at blow position of transponder surface under different crosswind velocity

由图16可知,不同横风风速下,应答器表面击打位置处的应力变化时程曲线线形相似,数值相近。为进一步提高数据的可靠性,减小随机性,在每种横风风速下进行10次击打,绘制的最大应力与横风风速关系散点图如图17所示。提取不同横风风速下,应答器表面受到的最大应力的平均值如表3所示。

表3 应答器表面的最大应力平均值Tab.3 Average maximum stress on transponder surface

图17 最大应力与横风风速的关系Fig.17 Relationship between maximum stress and cross wind speed

由表3可知,在行车速度为250 km/h的情况下,不同横风风速影响时,应答器表面所受到的最大应力平均值相差不大,范围为10.064～10.090 MPa,最大最小值之间仅差0.026 MPa。与无横风影响时的10.052 MPa相比,有横风影响时,风振耦合作用下应答器受力略有增加,但最大最小值间仅相差0.39%,说明横风作用对掉冰击打应答器的影响可忽略不计,应答器在冰雪击打下的受力主要还是随着高速运行的车辆掉落冰块产生的力所致。

3.3 冰块质量影响研究

冰块质量对于应答器的击打具有重要的影响,本节针对动车组车底部掉落冰块的大小进行探讨,根据现场实际情况,通过调整冰块的体积,得到3～30 kg按3 kg递增的不同冰块,探究在行车速度为350 km/h的列车空气动力效应和列车动荷载的风振耦合作用下,车厢底板结冰脱落击打应答器后,应答器所受最大应力随冰块质量的关系,如图18所示。

由图18可知,随着冰块质量的逐渐增大,应答器击打所受的最大应力呈现先迅速增大后缓慢增加的趋势。这说明随着冰块质量的增加,击打时的动能明显提升,应答器击打所受的最大应力随之提高,后续缓慢增加,这与增大了冰块与应答器击打时的接触面积密切相关。在冰块质量为30 kg时,最大应力为31.440 MPa,约为3 kg时的2.02倍。因此,采取相关措施降低覆冰体积从而减小击打的冰块质量,对于保护应答器的良好状态,维护列车安全运营具有重要的意义。采用对数函数可较好地表征出二者之间的相互关系,如式(4)所示,决定系数R2值为0.994,拟合效果较好。

(4)

4 结 论

针对高速铁路有砟轨道在冰雪条件下产生应答器被击打的问题,率先采用DEM-MFBD-CFD耦合分析方法进行了模型建立与影响因素分析。主要结论如下:

(1)基于严寒地区冰雪飞溅击打应答器问题,以CRH380B型动车组为例,创新性地建立了车厢底板结冰脱落击打应答器精细化仿真模型,实现了冰雪脱落击打应答器这一系列过程的准确模拟,并通过风洞试验,验证了模型的有效性。

(2)行车速度是应答器击打的重要影响因素。研究表明,应答器受到的最大应力与行车速度呈线性增长。为降低应答器受击打程度,在一些特殊地段、冰雪飞溅严重的情况下可以采取降速的措施,降速是降低冰雪飞溅击打危害的最有效的措施。

(3)风振耦合计算表明,横风作用对掉落冰雪击打应答器的影响可忽略不计,应答器在冰雪击打下的受力主要还是随着高速运行的车辆掉落冰块产生的力所致。冰块质量对于应答器击打呈现先迅速增加后缓慢增长的趋势。降低覆冰体积对于降低应答器击打具有重要意义,需提前做好融冰除雪工具设备、人员储备,规范使用除冰融雪工具,加强动车组除冰融雪作业。