基于POD分解的高速列车尾流动力学特性研究

2020-10-09郭迪龙张子健杨国伟

铁道学报 2020年9期

刘雯，郭迪龙，张子健，杨国伟

(1. 中国科学院力学研究所，北京 100190； 2. 中国科学院大学工程科学学院，北京 100049)

高速列车尾流区具有复杂的三维流动结构，由剪切层、分离与再附区、脱落涡系(尤其是一对大型的反对称旋转流向涡)等构成，而且这些脱落的涡系在远离列车的过程中会向下、向外运动，并且交替破裂[1-4]。尾流区的这些流动结构对高速列车产生了非常不利的影响，不仅增大了头尾车的压差阻力，也严重影响尾车的升力和侧向力，对乘客乘坐的舒适性和列车运行的稳定性、安全性等都是一大挑战[5]。列车风峰值出现在尾流场的近尾流区，会产生显著的诱导压力，这将给站台的乘客、设备和铁路沿线工人带来极大的安全隐患[2,6]。此外，对于更高速度的磁悬浮列车，其车体轻，气隙小，对升力和侧向力的要求非常高，一旦控制不好，就会增大车体与轨道的磨损。因此，研究列车尾流的非定常特性具有重要意义。

在列车尾流中存在着一对反对称旋转的流向涡，是许多类似流动所共有的特征，这是Ahmed在研究其建立的通用地面车辆中发现的[7]。近年来，高速列车的非定常尾流研究得到了国内外学者的高度重视。通过大量的研究[8-9]发现，IDDES(Improved Delayed Detached Eddy Simulation)方法在模拟压力系数分布、描述流场精细结构等方面均能取得较好的结果。Bell等[10]在1∶10的风洞试验中通过眼镜蛇探针(Cobra Probe)及表面压力测量方法，发现ICE3高速列车近尾流区存在周期性涡脱落，对应的无量纲频率(斯特劳哈数)St=0.2。Xia等[11]通过PIV(Particle Image Velocimetry)实验对CRH3高速列车尾流场中8个位置处的脉动速度进行了监测，得到其尾流场的主频St=0.18。Weise等[12]在风洞试验中采用油流和烟雾流动显示方法对不同的列车尾车形状进行了研究，发现尾流区会出现对称的分离泡或涡脱落的流动结构。Yao和Guo等[8]通过数值计算方法研究了尾流场在3个运行速度下的瞬时流场特性。然而，目前已有的诸多研究都是基于实验流动显示或仅对数值模拟得到的瞬时流场进行分析[8,11-12]，这是一种定性的分析手段，并不能很准确地对尾流场中的尾涡参数进行量化；或者是通过对某几个点的时间演化来定量分析尾流场的列车风和气动特性[10]，不能实现对整体流场的定量分析。此外，无论是实验还是数值模拟都还未能把尾涡结构彻底搞清楚，所以还需要更多深入的研究。

本征正交分解(Proper Orthogonal Decomposition, POD)方法是一种有效的流场降阶、简化方法，它将流场分解成一系列基函数和其随时间变化系数的线性叠加，并将这些基函数按照能量(即对流场的贡献度)进行排序。因此，通过捕捉高能量的流场模态，可以有效地提取流场的重要相干结构，从而更好地描述和剖析高速列车尾流场的演化特性。范晨麟等[13]和罗杰等[14]使用POD方法对多段翼流动进行了简单分析，发现POD方法对流场信息有着快速、有效地预测功能。Muld等[15]使用POD方法提取了ATM列车模型近尾流区的流场结构，成功识别到一对反向旋转的流向涡，求得该涡对脱落的无量纲频率St=0.085，并通过流场重构发现了该涡对在演化过程中的弯曲运动。因此，POD方法是进一步分析高速列车尾涡动力学问题的一个重要手段。

本文使用IDDES方法对CRH380A三编组列车尾流场的非定常特性进行了研究，分析运行速度200～450 km/h范围内列车尾涡的动力学特性，并通过POD方法对尾流场进行了模态分解，对模态结构、频率等流场特性进行了对比分析和讨论，从而得到高速列车尾流场在演化过程中的重要相干结构和主要规律。研究结果为以提高乘客乘坐的舒适性，列车运行稳定性、安全性为目标的尾车气动外形优化设计提供一定参考。

1 计算模型与计算方法

1.1 计算模型与计算域

本文以1∶1复杂外形的CRH380A三编组高速列车为计算模型，见图1，考虑车体风挡和转向架等重要结构，忽略受电弓、车顶空调导流罩等附属结构。以列车高度H为特征长度，车总长约22.3H。计算域尺寸为91.4H(长)×17.14H(宽)×12H(高)，来流距离头车鼻锥22.7H，出口距离尾车鼻锥46.4H，见图2。在求无量纲频率St时，以来流速度V作为特征速度，即St=f×H/V，其中f为频率。定义列车中心为x=0，y=0，车轮底端为z=0。定义流动方向为x轴正方向，竖直方向为z轴正方向，y轴正方向可由右手定则确定。

图1 CRH380A计算模型

图2 计算域示意

边界条件设置：由于列车运行速度450 km/h处于亚声速可压缩流动区，入口、出口以及四周均采用无反射的自由流边界条件。为精确地计算出列车表面的摩擦阻力及压力分布等，车体及车上各部件均采用无滑移壁面边界条件。通过设定列车静止、地面移动，从而考虑地面效应的影响。

1.2 网格拓扑

空间网格为正交切体六面体网格，边界层网格为三棱柱和金字塔网格。在壁面设置10层边界层网格，第一层边界层网格高度为1.71×10-4H。为了和六面体网格更好地衔接，并保证网格质量，设置边界层网格增长率为1.5。为保持网格的一致性，不同列车速度采用同一套网格进行计算。列车速度不同，相同的壁面第一层网格高度所对应的无量纲值y+会略有差异，但当前设置的壁面第一层网格高度均能使不同车速(200～450 km/h)所对应的y+值都很好地落在标准壁面函数最佳适用范围[16-17]之内(y+=30～150)。

y=0截面上和尾车表面的网格分布见图3。为了捕捉高速列车尾流场的更多细节并保证计算精度，在列车尾部区域进行了网格加密，保证体网格尺寸不超过1.79×10-2H。此外还对整车进行了3层网格加密，保证最小体网格尺寸不超过3.57×10-2H。网格总数量约为5 400万。

图3 网格分布示意

1.3 计算方法

列车近尾流区流场极其复杂，为了准确模拟流场中的细节，采用基于SSTk-ω湍流模型的IDDES方法求解可压缩N-S方程。目前，近壁湍流的处理有2种方法：直接求解到壁面和采用壁面函数。若近壁湍流直接求解到壁面，需使用y+< 1的网格。考虑到本文的计算模型为1∶1真实外形高速列车，其外部流动是复杂的三维流动，总网格量相当大。因此，为控制网格数量并参考众多学者的处理方法[8,18]，对近壁湍流的处理采用了标准壁面函数，并通过设置壁面第一层网格高度使y+值很好地落在最佳适用范围之内(y+=30～150)，从而保证了标准壁面函数应用于湍流边界层模拟时的准确性。

无黏项采用Weiss-Smith预处理的Roe差分方法进行离散，黏性项采用二阶中心差分方法进行离散。时间上采用二阶隐式方法进行积分。为了准确地模拟尾涡结构在每个周期内的发展情况，采用双时间步长法，设定Tref=H/V，全局时间步为Δt≈0.012 5Tref，其内迭代数为15，能够保证每个时间步内残差值至少下降2个数量级。以上计算方法通过STAR-CCM+实现。

2 本征正交分解(POD)

对于时刻i的瞬时流场，在nx个空间位置上进行采样，并按一定顺序将其排列成一个列向量vi，称为流场的一个快照。考察nt个时刻，即nt个快照，则可排列组成一个矩阵A，A=[v1,v2,…,vnt]。矩阵A代表了流场的空间-时间信息。在POD分解[15,19-20]中，流场被分解成一系列代表流场本征模态的基函数及其相应模态系数的线性叠加，即

( 1 )

式中：φ(k)(x)ak(t) 为一个流场模态;φ(k)(x) 为模态的基函数;ak(t)为模态时间系数，其具体表达式为

( 2 )

其中：λk、α(k)分别为矩阵C的第k个特征值和特征向量，且λk按降序排列。矩阵C为A的时间关联矩阵，即

C=ATA

( 3 )

Cα(k)=λkα(k)k=1,2,…,nt

( 4 )

λ1>λ2>…>λnt

( 5 )

3 POD方法验证

为了说明POD方法在提取流场相干结构方面的适用性和准确性，这里以比较经典的、已被大量研究的二维圆柱绕流为例，采用POD方法研究雷诺数Re=8 000时单圆柱后的尾迹流场特征。通过数值计算得到了非定常条件下圆柱升力系数随时间的变化曲线，见图4。大量研究表明[21-22]，当雷诺数Re处于亚临界区(300 ≤Re≤ 3×105)时，升力变化的无量纲频率St约为0.2。这里，通过对升力系数的时间序列做傅里叶分析(FFT)，得到的无量纲频率St=0.201，见图5。

图4 圆柱升力系数随时间的变化曲线

图5 频率特性示意

为了进行POD分析，在圆柱后侧10D×5D(圆柱直径D=0.01 m)的区域内均匀布置了33行、71列共2 343个监测点，提取流场在这些位置的2 000个时刻的速度分量u和v，得到2 000个流场快照。通过POD分解得到2 000个流场模态，并对各个模态时间系数做傅里叶变换，可以得到各个模态的频率特性。

由图5可见，对POD二阶模态的模态系数做傅里叶变换得到的无量纲频率St=0.201，与直接对升力系数做傅里叶变换所得结果完全一致。图6为通过POD分解后尾流区的前几阶模态图，与Zhang等[23]的试验研究成果一致。通过以上分析，证明了POD分解方法在提取流场相干结构方面的有效性。

图6 圆柱绕流POD模态

4 结果与讨论

4.1 瞬时流场尾涡特性

不同运行速度下高速列车的尾涡结构、涡脱落频率等流动特性不同，相应的尾车气动性能也不同。本节通过数值模拟，以运行速度200、300、400 km/h为例，给出不同速度下列车尾流场的压力分布、涡量分布和Q等值面图等，直观地表示出瞬时流场的尾涡结构形态及其发展规律，以便于进一步与POD分解的各阶流场模态进行对比，从而揭示高速列车尾涡动力学更深层次的规律。

不同运行速度下列车尾流场的压力云图和Q等值面图见图7。由图7可以看出：气流在尾锥处脱落出2个结构尺寸较大的涡，随后，这2个涡以螺旋状向下游发展，涡结构逐渐增大，并且在向下游发展的过程中与地面相互作用而破碎成许多强度不一的小涡。此外，整个瞬时流场的涡街结构类似于“树杈”形，在流向方向上越向下游，涡对间距越大，这是涡对之间的相互诱导作用所引起的。

图7 不同运行速度下列车尾流场的Q等值面图和瞬时压力云图(Q=150)

另外，由图7中还可看出，列车尾流中的涡街结构数量随着列车速度的增加并未发生改变，均约有4个完整的涡结构(图7中红圈所示)。通过简单的数学推导可得，相同范围内的涡结构数量与无量纲频率St成正比，由此可知3个速度下的St相差不大。一般情况，如果频率与速度呈线性关系，则St相等，以上结果正说明了涡脱落频率与运行速度之间呈线性关系。然而如前所述，这只是定性分析，不同速度下具体的涡脱落频率，将会在第4.4小节中通过POD方法进行进一步的定量分析。

高速列车复杂尾流场的形成与列车底部复杂结构以及地面效应密切相关，不同运行速度下沿高度方向不同截面的瞬时涡量图见图8。由图8可以看出：在地面附近，由头车第一个转向架发展出两股低速气流，该气流沿着流向迅速混合成一股较大的气流，并与风挡和其他转向架周围的气流相互作用，逐渐形成强度不一的旋涡。气流在流出尾锥之后，进一步与尾涡相互作用，形成更复杂的涡街结构。在离地面近1 m高的截面上，明显可见速度越高，尾涡的强度越大，其沿流向的影响范围也越大，尾流场越混乱。在离地面近2 m高的截面上，列车周围流场趋向于平稳，涡核强度减弱，列车底部结构和地面效应对尾流场的作用减弱。在离地面近3 m高的截面上，尾涡基本消失。通过比较不同高度下尾流场的涡量图，体现了列车尾涡的近地性。

图8 不同速度下不同高度截面处的瞬时涡量

4.2 POD模态分解

在数值模拟的基础上，本节采用POD分解方法对高速列车尾涡的时间演化规律做进一步研究。在尾流区选取距离尾车鼻尖1H～ 8H范围内均匀布置15个站位(间隔为0.5H)，每个站位上数据点为160×50个，提取流场在这些位置的速度、平均压力、脉动压力、涡量等进行POD分解。为了保证所提取模态的收敛性，需保证足够多的快照样本，以获得足够的周期数。参考Muld等[15]的参数设置和收敛性验证方法，设定采样时间步Δtpod≈0.25Tref，总样本时间TA≥ 250Tref，该样本能够保证至少涵盖45个周期，且所得到的前四阶模态的结果是收敛的。

不同运行速度下列车尾流场的3个速度分量u、v、w的前四阶POD模态见图9(速度等值面，蓝色为正值，红色为负值)。一阶模态是平均流，这是以平均流作为主导的流动类型的一个POD特性。同时，一阶模态也是最主要的模态，在POD分解的各阶模态中能量最高，比其他模态高出至少3个数量级，见图10。通过分析速度分量v和w的等值面可以看出，一阶模态代表的是高速列车尾流中的2个对称的、反向旋转的流向涡。在w的一阶模态等值面图中，中心周围有强烈的向下流动，区域两侧有强烈的向上流动；在v的一阶模态等值面图中，左半区域的上侧为向右流动，下侧为向左流动，右半区域与此相反。综合来看，左半区域中上面向右流动，在区域中心转为向下流动，下面向左流动，在区域左侧再转为向上流动，因此形成一个顺时针旋转的涡旋；同理，右半区域形成一个逆时针旋转的涡旋，见图11(a)。通过观察w的一阶模态中流动结构的尺寸和位置发现，沿着流向两涡间距逐渐增大，与4.1节瞬时流场的结论一致。

图9 不同运行速度下列车尾流场的前四阶POD模态

图10 各阶POD模态的能量占比示意

二阶及以上模态代表一系列非定常振荡。二阶模态展示出的是涡对的交替生成与脱落这一流动状态，三阶模态与二阶模态是成对存在的，同样代表着涡对的交替生成与脱落，但所对应的旋涡的旋转方向与其相反，即存在着相位差，见图11(b)、图11(c)，Muld等[15]认为这种关系对应于流动结构中的对流。比较四阶模态与二阶模态发现：四阶模态与二阶模态类似，但四阶模态中比较清晰的结构数量约有4个，而二阶模态中仅为2个，且四阶模态的空间结构更小，说明四阶模态代表的非定常振荡频率更高，约是二阶模态的2倍关系，见图11(d)。

图11 模态几何示意

比较不同运行速度下的POD各阶模态发现：不同速度下同一模态的整体空间结构是相似的，结构数量和尺寸大小也几乎一致。以四阶模态为例，3个速度下在相同的范围内均约有4个完整的流动结构，这与4.1节中瞬时流场的涡结构描述是一致的，证明涡脱落频率与速度呈线性关系，进一步的定量分析将在4.4节中讨论。这种线性关系也理应反映在二、三阶模态的流动结构中，但由于二、三阶模态数量相对四阶模态少，不易简单地从模态流场结构图中识别出具体差别，有待进一步的定量分析。

4.3 POD对流动频率的解析能力

在对流场的动力学特性进行分析时，尤其是流场的频率特性，不管在实验中还是数值计算中，往往很难实现全局的定量分析。比如，实验或数值计算中往往通过对局部几个测点进行测量(比如压力传感器，热线风速仪，或其他监测手段)并进行傅里叶分析来得到流场的动力学特性，虽然是定量的分析，但却是局部的，以此局部量来代表流场整体的流动特性往往需要对流场更深层次的认识和慎重处理；又或者常常通过流动显示技术(比如烟迹法、PIV等)、瞬时流场分析(比如压力云图、涡量云图)等手段，但这却仅仅局限在定性分析上，不能实现对流场动力学特性的量化分析。

在对尾涡的非定常特性进行数值模拟分析时，传统的方法往往是在列车表面及尾流区布置压力测点，通过对测压点的压力随时间的变化做傅里叶变换，从而得到尾流区的周期性变化特性，即尾涡脱落的频率等。为了说明这种方法的局限性，以运行速度300 km/h为例，分别在尾车表面布置3个压力测点(①、②、③号)，尾流区布置3个压力测点(④、⑤、⑥号)。监测点的布设位置见图12。以尾车鼻尖在地面的映射点为原点，①～⑥号测点的位置(X/H，Y/H，Z/H)分别为(-0.07，0.06，0.38)、(-1，-0.33，0.54)、(-0.66，0.3，0.46)、(1，-1，0.3)、(3.2，0.5，0.54)和(1，0，0.3)。

图12 压力测点示意

通过监测1～6号测点在20 000个物理时间步(250Tref)内的压力变化，并对其进行傅里叶变换，得到每个测点所对应的功率谱密度图，见图13。由图13可以看出，列车表面及尾流测点压力的频谱结果频带较宽，且杂乱，没有明显的峰值，这给尾流场的频率分析带来较大的困难，不利于识别尾涡脱落频率等，也不利于进一步对振动与噪声进行相关分析。此外，通过单一测点得到的结果是局部的，仅与当地流动的非定常特性有关，不同测点之间的频谱差异较大，很难据此给出整体流场的频率特性，而往往整体流场的振荡特性才真正与尾车振动密切相关。

图13 测点①～测点⑥的傅里叶变换结果

POD方法可以很好地弥补上述问题的不足，其通过分析整个流场的时间相关信息，对我们关注的尾流场区域进行整体的降阶分解，并根据各个模态的能量占比来遴选出各阶主模态，并得到各阶模态所对应的频率。这样得到的各阶模态之间是近似解耦的，每个模态具有较明显的主频，大大提高了频率分析的准确性。POD二阶模态所对应的频谱，见图14。

图14 POD二阶模态的频率特性

4.4 POD模态频率特性

高速列车尾流的非定常特性是由尾车上交替的涡脱落产生的，涡脱落的交替性导致尾流具有一定的周期性。通过POD方法，高速列车尾流场被分解成一系列随时间变化的模态，每一个POD模态(一阶模态除外)代表了尾流的一种周期性运动。图15(蓝色为正值，红色为负值)通过选取速度分量v在一个周期内的12张结构差异比较明显且具有代表意义的模态快照，给出300 km/h速度下高速列车尾流场POD二阶模态在一个完整周期内的演化过程。

图15 300 km/h速度下v的POD二阶模态的一个运动周期

通过图15可以看出，快照223至234中，尾车附近是一个逆时针的涡旋，且强度先增大后减小；快照237至247中，尾车附近变成了一个顺时针的涡旋，强度同样先增大后减小。显然，随着时间的推移，这样一对方向相反的大型流向涡在列车尾部交替地生成与脱落。二阶模态的这一流场特征代表了列车尾流的主振荡模式。通过快照所显示的流场结构可以看出，二阶模态是随时间周期性变化的，且周期不是严格均匀的。统计40个周期共有约1 140张快照，一个周期平均约28.5张快照。POD采样时间步长为10 ms，则平均周期约为285 ms，对应的物理频率为3.51 Hz。通过特征长度H和特征速度V作无量纲化，得到该速度下二阶模态的无量纲频率St=0.147。

如前所述，POD方法在分析频率特性方面具有优越性。每个POD模态的模态系数ak(t)是随时间变化的，这启发我们可以通过对模态系数做傅里叶变换来分析每个POD模态的频率特性。以运行速度300 km/h为例，图16给出了二阶、三阶、四阶模态的FFT结果。该速度下二阶、三阶、四阶模态具有明显的主频，其对应的St分别为0.148、0.148、0.268。