基于延迟脱体涡算法高速列车通过隧道时的绕流特性

杨广雪，刘志明，刘秋泽，董 磊，周君峰

（1.北京交通大学 机械与电子控制工程学院,北京 100044；2.中车长春轨道客车股份有限公司 转向架开发部,吉林 长春130062；3.中车青岛四方机车车辆股份有限公司 技术中心,山东 青岛266111）

高速列车通过隧道时：由于气动效应和地面效应的共同影响,绕流更加复杂,对列车及周围环境产生更严重的影响；列车运行阻力会明显增大,甚至可能高出明线运行的4倍以上,对列车的安全运行提出了更高的要求［1］；在隧道和列车之间会产生强烈的压力波动,诱发列车内压力波动加剧,降低车辆的乘坐舒适性［2-3］；由于列车运行轨迹偏离隧道中心,车体两侧的压力波动差增大,使得列车横向振动加剧,从而降低列车的运行稳定性［3］；列车在通过隧道出口、入口时车身结构受到的压力幅值可分别高达4 和7 kPa 左右,远远超出明线运行时的压力幅值,这会对列车相关结构的承载强度造成更大的影响［4］；列车在进入隧道时的挤压效应会在隧道的出口部位产生压力波动,对周围的环境产生噪声污染和振动影响［5］。基于这些原因,进行高速列车过隧道绕流特性研究具有重要的理论意义和工程应用价值,但国内外针对高速列车通过隧道时的绕流特性研究仍缺乏较为全面的分析。

目前,高速列车绕流的模拟方法主要包括：直接数值模拟（DNS）、大涡模拟（LES）、雷诺平均方法（RANS）和脱体涡模拟（DES）［7-10］。DNS 方法通过直接求解流体运动的N-S 方程,计算精度最高且结果最详细,但对网格分辨率和计算资源要求很高,就目前的计算机水平而言,直接将DNS 方法应用于高速列车的绕流模拟仍不现实,除非是在极其简化条件下才是可行的；周瑜［11］通过求解减少一个维度的二维非定常黏性可压N-S方程,对极其简化的高速列车模型通过隧道时的流场进行了数值模拟,得到了不同时刻系统二维流场云图。LES 方法通过先对大尺度涡和小尺度涡过滤,然后分别对小尺度涡用模型求解和对大尺度涡直接数值求解,该方法求解精度较高,消耗资源比DNS 方法少；鲍龙［12］采用LES 方法对不含转向架的ICE2和CRH380A 型高速列车在复线隧道内行驶时的外部绕流进行了数值模拟,得到了列车绕流的瞬时速度及压力云图。RANS方法通过将N-S方程进行时间平均,将非定常的湍流问题转化为一个定常的问题进行研究,该方法可以降低时间和空间的分辨率,减少对计算系统的需求并能获得符合工程要求的计算结果；赵晶［13］基于雷诺平均N-S 方程及两方程紊流模型,对CRH2型动车组过隧道时的空气动力学问题进行了研究。DES 方法在湍流附面层内采用RANS 方法而在其它区域采用LES方法,既节省计算成本而且计算结果也比较令人满意,但在建立模型时处理不当会引起模型的应力损耗。为此,Spalart 等人［14］提出了延迟脱体涡算法（Delayed Detached Eddy Simulation,DDES）,即在DES 计算中加入参数延迟控制,从而有效解决DES 方法容易引起应力损耗的问题,逐渐在列车外部绕流的数值模拟中得到广泛应用。

因此,本文结合目前国内外研究现状及存在的问题,基于延迟脱体涡算法和滑移网格技术,建立CRH380A 型高速列车含有转向架的三维可压缩瞬态仿真模型,模拟研究高速列车通过隧道时的气动力、速度场和表面压力3 大绕流特性的变化规律,以期为高速列车相关结构的强度评估和设计提供较为可靠的数据参考。

1 仿真模型的建立

1.1 延迟脱体涡算法

本文采用的DDES算法是对Spalart 等［14］于1997年提出的脱体涡算法的改进,即通过在边界层内进行参数延迟控制,从而最大程度地降低DES算法引起的应力损耗。

DDES 算法的基础是Menter 于1994年提出的剪切应力输运（SST）k-ω湍流模型方程［15］,即

式中：ρ为流体密度；k为湍流运动的湍动能；t为时间；Ui为沿i方向的平均流速分量；xi为i方向的位移；Pk为湍流应力；μ为流体动力黏度系数；μt为涡黏系数；ω为湍流比耗散率；F1为混合函数；β*,σk,γ,β,σw和σw2均为模型常数,不同计算条件对应不同的取值。

混合函数F1的定义式为

其中,

式中：y为近壁面最短的距离；Dt为交叉扩散项。

延迟脱体涡模拟算法的核心是将k方程式（1）中的耗散项（βρkω）修改为

其中,

式中：M为网格中心到相邻单元中心距离中最大的一个；CD,CD,k-ω,CD,k-ε均为模型常数；F2为另一混合函数；z为z方向的位移。

混合函数F2的定义式为

式中：ν为流体运动黏度系数。

在延迟脱体涡算法中：若FD= 1,则采用RANS模型计算；若FDES＞1,则采用LES模型求解。如果令F2= 0,则为原始的脱体涡算法。已有文献表明［16］,在原始的脱体涡算法中,模型系数CD采用0.65 是最优值；而在延迟脱体涡模拟算法中,模型系数CD,k-ω=0.78,CD,k-ε=0.61。

1.2 滑移网格

基于滑移网格技术实现数值模拟中列车与隧道之间的相对运动。如图1所示,将计算域分为区域1 和区域2 两部分,区域1 为列车所在空间网格,区域2 为隧道内列车以外的空间网格。2 个区域独立划分网格,相邻的网格则形成1个交界面；计算过程中,区域1 和区域2 的流场信息通过这一交界面进行交换。以区域1中的单元Ⅰ为例,单元Ⅰ的信息通过交界面a-f和f-b传递到区域2 中的单元Ⅴ和Ⅵ；其他单元信息交换与此类似,从而实现流场的瞬态数值模拟。

图1 滑移网格示意图

1.3 仿真模型及定解条件

列车模型基于高速动车组CRH380A的1/25缩比简化模型,包括1 节头车和1 节尾车（含司机室）,头车的2 个转向架采用简化转向架模型,尾车转向架采用相对真实的转向架模型,采用如图2所示。列车模型长度ltr=1.58 m。

隧道模型参考京沪高速铁路大断面隧道［17］,半径为6.87 m,横断面面积为100 m2,双轨线间距为5 m ,列车断面积的阻塞比为0.12。隧道长度ltu是列车长度的5 倍,即ltu=5ltr,从而尽可能准确地捕捉列车实际受到的气动特性。

整个计算域的几何模型及边界条件如图3所示,图中H代表列车高度,作为特征长度。计算域包括隧道域、隧道出入口的空气域和列车周围的移动域；隧道壁面、列车表面边界均设定为无滑移壁面,其余边界条件均设定为压力远场。流场信息则通过将隧道空气域与列车移动域的边界设定为交界边界来传递。整个计算域的离散方式采用六面体网格,并在车体流场处进行网格加密,同时在转向架和车身表面处的位置采用边界层网格,如图4所示；隧道域的网格尺寸（方形网格边长）为6.8 mm,列车移动域的网格尺寸为3.4 mm,设置边界层内的网格增长率（网格尺寸的变化率）小于1.3,边界层外的网格增长率小于等于2。

图2 基于CRH380A的简化列车模型及其转向架模型

图3 几何模型及边界条件

图4 纵截面的计算网格

根据相关文献［18-19］,当Ma2ltu/ltr≪1（Ma为马赫数,Ma=Ui/c,c为声速）时,在模拟计算列车的绕流时必须考虑列车进出隧道时引起的空气可压缩性。基于此,建立基于延迟脱体涡算法的三维可压缩非稳态仿真模型。根据Sutherland 公式,理想气体的动力黏度系数μ计算公式为

式中：S为Sutherland常数,本文取111 K；T0为参考温度,其值为293 K；μ0为参考温度下的动力黏度系数,其值为1.79×10-5Pa·s。

计算流场的初始化：设置压力远场U∞=0；列车以25 m·s-1的速度运动,以此模拟列车进出隧道的过程；计算初始,隧道入口与列车相距约是列车长度的1.5 倍,超过文献［12］中给定的隧道入口与列车相距20～30 m,将有利于促进湍流的充分发展；时间步长设定为0.000 1 s。依据列车25 m·s-1的运行速度,列车与隧道刚接触的时刻为t1*=t×Vtrltr= 1.38 s,列车刚全部进入隧道的时间为t2*= 2.38 s,列车尾部与隧道入口相距列车长度的2倍在t3*=4.38 s。

2 气动力

根据列车阻力主要来源于列车的头尾压差阻力,为了便于分析,将列车划分为3个部分：头车Headbody、尾车Tailbody 和车底Others,其中车底包括列车底部所有结构,转向架及设备舱、设备舱底板、排障器底板）。

2.1 气动阻力

图5 列车进入隧道过程中气动阻力随时间的变化曲线

列车进入隧道过程中气动阻力随时间的变化曲线如图5所示。图中：tH为头车到达隧道口的时间点；tT’为车尾刚驶入隧道的时间点；tT为车尾完全进入隧道的时间点；Fd-Head,Fd-Tail,Fd-Others,Fd分别为头车、尾车、车底和整车的气动阻力。由图5可看出：列车头部刚驶入隧道时,头车的气动阻力突然增加,整车的气动阻力也随之骤升；当车头完全进入隧道时,头车的气动阻力也达到最大；当尾车刚驶入隧道时,头车的气动阻力逐渐降低,整车的气动阻力也随之减少,尾车的气动阻力逐渐增大；当尾车全部刚刚进入隧道时,尾车的气动阻力骤然再次增加；当列车全部驶入隧道内时,尾车与头车的气动阻力呈现出相反相位的气动波。

2.2 气动升力

列车进入隧道过程中气动升力随时间的变化曲线如图6所示,图中,Fr-Head,Fr-Tail,Fr-Others,Fr分别为头车、尾车、车底和整车的气动升力。由图6看出：在头车刚进入隧道时,列车3个部分的气动升力均无明显变化；随着列车的不断进入,大约在尾车刚驶入隧道时,列车3个部分的气动升力均明显升高；而当尾车全部刚刚进入隧道时,列车3个部分的气动升力无较大变化并且呈周期性变化；当列车全部驶入隧道内时,头车和尾车分的气动升力均为正,这有利于减小轮轨力,而车底的气动升力为负,将列车压向轨道。

图6 列车进入隧道过程中气动升力随时间变化曲线

2.3 气动侧向力

列车进入隧道的过程中气动侧向力随时间的变化曲线如图7所示,图中Fs-Head,Fs-Tail,Fs-Others,Fs分别为头车、尾车分、车底和整车的气动升力。由图7看出：车底的气动侧向力无明显变化,保持在0 上下小幅波动；头车和尾车的气动侧向力变化趋势基本一致,头车的波动幅度和数值均比尾车大；在尾车将要进入隧道时,头车和尾车的气动侧向力突然增大,且均为正值,表明气动侧向力使得列车被压向靠近隧道的一侧。

图7 列车进入隧道过程中气动侧向力随时间变化曲线

3 速度场

为了观察列车通过隧道整个过程中车下流场的变化规律,在列车下方距离轨面z=-0.025 m 的高度上选取4 个横断面,其位置分别为距车头x1=7.49 m,x2=16.49 m,x3=26.00 m 和x4=36.49 m,在每个横断面上沿车宽方向的位置设置3～7 个测点,测点分布在y=±2.25 m 之内,共21 个测点,编号为Lat 1—Lat 19,Lat 40,Lat 41,Lat 42,详细如图8所示,通过设定,使得这些测点相对于列车的位置始终不变。由于在列车高速运动前提下,测点的速度变化比较剧烈,所以采用移动平均方法对测点的速度数据进行处理,处理时将时期数（用于求平均数的数据量）设定为100。

图8 车下测点位置（单位：m）

3.1 纵向速度分量

在列车进入隧道过程中不同测点的无量纲速度纵向分量vx随时间的变化曲线如图9所示,由图可得如下结论。

（1）在y=±2.25 m 线上,列车两侧对称测点的速度纵向分量变化趋势一致；头车测点Lat 1 和Lat 6 的速度纵向分量,在头车进入隧道时急剧增加,之后随着头车不断进入隧道,其值基本无变化；另外3个断面上测点的速度纵向分量,在头车到达隧道时开始减小并且流动方向逐渐反向,当尾车进入隧道口时基本保持不变,在尾车进入隧道口时隧道与车身侧墙之间形成回流区。

（2）在y=±1.75 m 线上,各个测点的速度纵向分量变化趋势基本一致；当头车进入隧道时各测点的速度纵向分量逐渐减小但波动幅度增大；当尾车进入隧道时速度纵向分量近似在同一水平波动；处于对称位置的lat 8 测点,受隧道影响显著,其速度纵向分量变化相比lat 11测点更明显。

（3）在y=±1.00 m 线上,除了Lat 15 之外其它测点的速度纵向分量均有所降低,所有测点的速度流向未发生变化；位于靠近设备舱位置的测点Lat 9 和Lat 10,受隧道效应影响显著,速度纵向分量降低更明显,而其他测点速度纵向分量波动较为明显。

3.2 横向速度分量

在列车进入隧道过程中不同测点的无量纲速度横向分量vy随时间的变化曲线如图10所示,由图可得如下结论。

在y=2.25 m 线上,Lat 6 测点的横向速度分量在列车驶入隧道后明显减少,而Lat 13测点却有一定程度升高；而在y=-2.25 m 线上,各个测点的速度基本无变化。

在y=1.75 m 线上,测点的横向速度分量略微高于y=-1.75 m 线上的测点；在列车进入隧道过程中,隧道效应对各个测点没有造成明显的影响。

4 列车表面压力分布

高速列车表面压力分布情况是高速列车通过隧道时受到的气动效应的直接体现。为此,在列车靠近隧道壁的一侧,选取列车侧面、车身裙板和列车底部分别布置测点,分析高速列车通过隧道时其表面的压力分布。

图9 列车进入隧道过程中不同测点无量纲速度纵向分量vx随时间变化曲线

4.1 列车侧面压力分布

列车侧面上,在高度z=1.00 m、横向位置分别为x=1.7,7.5,11.5,15.0,23.5,36.5 m 处布置6个测点,如图11所示。

在列车进入隧道过程中列车侧面各测点压力随时间的变化曲线如图12所示,因S2与S3测点的变化规律十分接近,故未绘制S3 测点的压力变化曲线。由图可知：当头车到达隧道口时,5 个测点的压力都突然增大,由负压上升为正压,其中头车的2个测点S0和S1的增幅最大,其它测点增幅较小；当头车进入隧道之后,各个测点的压力开始下降,直到尾车完全进入隧道时才停止下降；此后,压力基本在负压范围内呈现周期性波动,处于尾流区的S5测点具有更强烈的波动特征。

4.2 列车裙板表面压力分布

图10 列车进入隧道过程中不同测点无量纲速度横向分量vy随时间变化曲线.

图11 车身侧面测点布置

图12 列车进入隧道过程中不同测点的压力随时间变化曲线

在靠近隧道壁一侧的车身裙板上,在高度z=0.63 m、横向位置分别为x=1.7,7.5,11.5,15.0,23.5,36.5 m 处布置8 个测点,如图13所示。

在列车进入隧道过程中裙板各测点压力（含头车4 个测点的压力标准差以及车身及尾车4 个测点的压力标准差）随时间的变化曲线如图14所示,由图可知：头车测点的负压呈现先减小后增加的趋势,当尾车完全进入隧道时,压力停止增加,此后各个测点的压力维持在较高幅值波动变化；位于车身和尾车的4 个测点负压在列车到达隧道口时开始大幅增加,在列车尾部到达隧道口时停止大幅增加,此后随着列车的前进,各个测点的压力维持在较高幅值波动变化；与列车未驶入隧道相比,列车驶入隧道之后头车4 个测点的压力标准差明显减小。

4.3 列车底部表面压力分布

在列车底部,即y=0 线上,选取7 个断面x=1.7,7.5,11.5,15.0,23.5,36.5 m 处,每个断面上布置1～3 个测点,共15 个测点,如图15所示。

图13 裙板测点位置

图14 列车进入隧道过程中裙板各测点压力随时间变化曲线

在列车进入隧道过程中列车底部各测点压力随时间的变化曲线如图16所示（为了便于分析,本文仅给出e1、e2、e4、e10、e13 和e15 这6 个测点的压力变化）。由图可知：这6 个测点的变化趋势具有相似性,列车刚进入隧道时测点的负压值出现略微的下降,之后随着列车的不断进入测点的负压继续下降,大约在尾车开始进入隧道时压力值迅速增加并在尾车完全进入隧道时停止增加,当列车完全进入隧道之后测点保持较高的负压水平上下波动；当列车完全进入隧道之后,6 个测点中,e15的压力水平明显低于其他测点,而e4的压力最高且呈现出同时存在低频振动和高频振动的复杂波动。

图15 车下测点位置

图16 列车进入隧道过程中列车底部各测点压力随时间变化曲线

5 结 论

（1）延迟脱体涡算法能较好地捕捉列车通过隧道时的气动特性；列车头部刚与隧道接触时,气动阻力迅速升高并在车头完全进入隧道时达到最大值；当尾车全部刚刚进入隧道时气动升力和侧向力骤然增加；当列车全部驶入隧道内后,气动力的波动幅值均明显升高；这均会对列车的平稳性和舒适性产生一定影响。

（2）在列车头部刚与隧道接触时,列车下方两侧的速度纵向分量会急剧增加,位于靠近设备舱位置的速度纵向分量会显著降低；当尾车刚与隧道接触时,隧道内壁与列车侧面之间的流场会出现回流区；而这会对隧道中的人身安全及设施安全产生不良效果。

（3）列车通过隧道时,列车侧面压力整体上呈现先增大后减小、最后维持周期性波动的趋势,处于尾流区的车尾部位具有更强烈的波动特征；列车裙板和车底的表面压力整体上均呈现出先减小后增加、最后维持在较高幅值波动的趋势。这会对列车相关结构的疲劳强度产生不利影响。