近年来，我国铁路发展迅速，随着路网持续的优化，高速铁路在高寒地区持续布局建设与运营使用是大势所趋. 寒区冬季低温多雪的气候，给高速铁路的发展带来新的挑战^[1]. 高速列车在积雪路线行驶时，沿线积雪受列车运行产生的近地表紊流影响，将克服颗粒自身的黏结与惯性影响，随紊流进入车底，在转向架部位形成覆雪^[2].

为了解决上述问题，学者们通过研究与实践形成了主被动结合的高速列车覆雪处理方案. 其中，被动除雪指应用水淋、热烘、敲击等技术^[3]实现车底结构表面融雪. 主动防雪指通过改变车底结构外形，如增设裙板，对列车与近地面空间的流场进行人为干预，降低雪颗粒随气流进入的频率，以减少车底的覆雪，逐渐成为国内外该领域的研究热点.

Trenker等^[4]基于标准剪切应力输运（SST）k-ω湍流模型求解列车周围流场，结合离散相模型（discrete phase model, DPM）对雪颗粒轨迹进行跟踪. 通过自定义模型壁面覆雪判别准则，包括壁面阈值摩擦速度与雪颗粒捕捉角度两方面，实现对模型壁面覆雪的可视化. Paradot等^[5-6]基于在Jules Verne风洞上开展的几何体覆雪测试，采用Trenker等^[4]的模拟方法与准则参数，开展模拟验证. 蔡路等^[7-8]利用DPM模型配合自定义覆雪判别准则，对高速列车转向架区域的雪颗粒运动与沉积结果开展模拟，指明了转向架各部件的平均覆雪体量. 张乐等^[9]探究该方法中的雪颗粒参数（密度与粒径）对覆雪结果的影响. Wang等^[10-13]基于DPM模型开展转向架、列车整车的覆雪模拟，讨论流场卷携雪颗粒进入车底的运动特征，提出多种防雪优化设计. 韩运动等^[14]对转向架流场状态进行实测，配合数值模拟，探究转向架覆雪的成因与气流的关系. 丁叁叁等^[15-17]对高速列车转向架开展流场模拟分析，从流场特征角度讨论转向架区域覆雪的成因及潜在分布.

总体来说，采用DPM模型对高速列车转向架区域覆雪的模拟尚不成熟，缺乏对壁面覆雪判别准则的细致探究. 本文基于Paradot等^[5-6]在JV风洞里开展的模型覆雪测试结果，对Trenker壁面覆雪判别准则^[4]进行讨论及合理改进. 利用DPM模型配合改进后的覆雪判别准则参数，对我国某型号高速列车的转向架模型开展覆雪模拟，明确了转向架易出现不利覆雪的区域，可以为高速列车的运行防雪优化提供技术支持.

高速列车转向架区域的覆雪属于典型风雪耦合作用的结果. 对于风雪耦合作用的模拟研究，针对雪相的模拟处理，可以划分为欧拉方法和拉格朗日方法. 对于前者，雪相被处理为与流体相（空气）相互贯穿的连续介质，通过引入相体积分数的方式区分控制体内的气固相，求解各相守恒方程. 对于后者，流体（空气）相为连续相，雪相被处理为离散颗粒系统，通过相间体积分数参数实现气固相的耦合，计算相间质量、动量的交换. 其中，以DPM离散相模型^[4-14]为代表的拉格朗日方法被广泛地应用于对高速列车转向架的覆雪模拟研究.

对于列车外部的流场模拟，采用不可压缩流体模型，流体连续性与动量守恒方程为

式中： $\;{\rho _{\rm{f}}}$、 ${\boldsymbol{u}}$与 $p$分别为流体的密度、速度与压力； $t$为时间； ${\boldsymbol{g}}$为重力加速度； ${\alpha _{{\rm{f}}} }$为流体体积分数， ${\alpha _{\rm{f}}} = 1 - {{{V_{\rm{p}}}} \mathord{\left/ {\vphantom {{{V_p}} {{V_{c} }}}} \right. } {{V_{{\rm{c}}} }}}$，其中 ${V_{\rm{c}}}$、 ${V_{\rm{p}}}$分别为控制体单元总体积及颗粒体积； ${{\boldsymbol{F}}_{{{\rm{f}}} - {{\rm{p}}} }}$为空气相与雪相间的作用力，主要为颗粒拖曳力 ${{\boldsymbol{F}}_{{\rm{drag}}}}$； ${\boldsymbol{\tau }}$为流体黏滞应力张量，

其中 ${\boldsymbol{I}}$为单位张量， ${\mu _{{\rm{f}}} }$为流体黏度.

流体黏度包括运动黏度与湍流黏度，可以在选用的湍流模型中应用，如RANS雷诺平均模型、LES大涡模型与DNS直接模拟模型等.

对于雪相，雪颗粒被视为球体，在不考虑颗粒间碰撞的前提下，力平衡方程为

式中： ${m_{\rm{p}}}$、 ${{\boldsymbol{u}}_{\rm{p}}}$与 $\;{\rho _{\rm{p}}}$为雪颗粒质量、速度与密度； ${\boldsymbol{F}}$为附加外场力，如电磁场力，此处可以忽略； $ {{\boldsymbol{F}}_{{\rm{drag}}}} $为雪颗粒在流场中所受的拖曳力，

式中： ${d_{\rm{p}}}$为雪颗粒粒径； ${C_{\rm{D}}}$为颗粒拖曳力系数，可以按光滑球模型^[18]，由颗粒雷诺数 ${{R} _{{\rm{e,p}}}}$所决定的一组经验常数 ${a_1}$~ ${a_3}$来确定：

采用传统SIMPLE方法求解流体（空气）相的控制方程，配合显式时间积分法，可以实现对颗粒在流场中运动方程的求解.

雪颗粒跟随流场输运与壁面发生碰撞，颗粒是否会被黏附在壁面，依赖于碰撞特性及壁面流场作用. 从能量的观点出发，雪颗粒与壁面的碰撞属于非弹性碰撞，碰撞后的动能损失程度决定了黏附或反弹行为. 雪颗粒动能的损失（消耗或转化为黏结内能）与雪自身的晶体类型、含湿特性、晶体间黏性等诸多复杂因素有关，无法用传统碰撞恢复系数或击溅函数来描述. 采用阈值摩擦速度 ${u _{*{\rm{t}}}}$来宏观描述壁面雪层的特征^[19]，如黏结、烧结、压实等. 依据 ${u _{*{\rm{t}}}}$与壁面切应力 $\tau $（摩擦速度 $ {\mu _*} $）间的关系，判别雪颗粒是否能稳定或持续黏附在壁面，即 $ {u _*} $> ${u _{*{\rm{t}}}}$时，壁面雪颗粒会被流体侵蚀，无法形成覆雪.

Trenker等^[4]提出壁面覆雪模拟判别准则：仅有同时满足 ${{{u}} _*} < {{{u}} _{*{\rm{t}}}}$与雪颗粒-壁面碰撞角 $\alpha < $雪颗粒阈值捕捉角 $ {\alpha _{\rm{t}}} $时，雪颗粒能够被壁面捕捉，以实现壁面的覆雪. 否则，雪颗粒将在壁面发生反弹与被侵蚀的行为，如图1所示. 其中，对 ${{{u}} _{*{\rm{t}}}}$的数值估算采纳了Nishimura等^[20]风洞边界层的风雪联合试验标定结果，取值为0.21 m/s. 该试验中雪颗粒经历了长期（达1 a）的冷库储存与水汽处理，晶体颗粒表面能大大降低，与真雪物性的差异较大. 蔡路等^[7-8]依据Kind^[21]的建议，对 ${{{u}} _{*{\rm{t}}}}$直接取实测估算极值1.0 m/s.

相比而言， ${u _{*{\rm{t}}}}$在不同的研究中取值差异高达5倍. 此外，类似参数的取值不确定性在 ${\alpha _{\rm{t}}}$中同样存在，文献[4~6]中45°取值的合理性亟待讨论.

壁面覆雪判别通过在模拟中设置自定义条件捕捉壁面边界实现. 利用ANSYS Fluent 17.0软件中的自定义存储器（user defined memory UDM），除参数 ${u _*}$与 $\alpha $外，对雪颗粒-壁面碰撞发生单元上的碰撞数 $n$、捕捉数 ${n_{\rm{t}}}$与粒径等参数进行记录. 模拟流程和UDM的嵌入应用与文献[8]一致.

为了保障计算结果的收敛性，先计算稳态流场至迭代收敛，后引入DPM离散相模型，包括颗粒入射、利用UDM记录颗粒-壁面碰撞特性与颗粒捕捉计数更新等方面.

为了系统地讨论基于DPM离散相模型的壁面覆雪模拟方法的可行性及准则参数的合理取值，选用Paradot等^[5-6]开展的楔体覆雪试验作为原型，开展同边界覆雪模拟方法的验证.

该原型试验在南特Jules Verne气象风洞中开展. 风洞为回流型，试验段尺寸为21 m×10 m ×7 m，最大试验风速为40 m/s，可以实现最低−25 ℃的低温环境，以开展风雪联合试验. 试验雪颗粒由造雪枪提供，通过调节喷射压力，可以实现对雪颗粒物性（密度、含水率）的控制. 具体的雪颗粒晶型为谷粒雪，如图2(a)所示，颗粒平均粒径为250 μm，对应的雪质量浓度为320~460 kg/m³.

如图2(c)所示，造雪枪被放置在风洞上游（x = 3 m），通过逆向喷射的方式将雪颗粒引入流场，颗粒受来流的影响，被吹至下游（x = 12.5 m）的测试平台，实现对平台上固定模型的覆雪. 试验前，通过在测试平台前放置雪收集盒，实现对雪质量浓度的标定. 结果显示，在风速为20 m/s与测试时长为10 min的试验条件下，测试平台前y向单位长度的雪质量浓度平均值为0.011 kg/m³.

测试模型被固定在测试平台顶板底面，是一组斜面坡度均为36° 的楔体（记为A与B），顶面为边长为0.2 m的正方形. 其中，模型A顶面垂直于来流方向（x向），模型B顶面的法向与来流方向成18°，考察不同斜面角度对覆雪结果的影响.

以楔体模型A为例，受模型布置与颗粒引入方式的影响，楔体迎风壁面均存在覆雪情况，且分布较均匀. 相比而言，楔体顶面覆雪累积量（或雪深）约为斜面的1.46倍，如图3所示.

基于原型，建立同尺寸、同边界模拟分析模型，包括一组楔体模型和用以模型固定的测试平台，如图4所示. 计算域尺寸为20H×10H×7H，其中H为固定架宽度（1 m）. 模型采用结构化网格，最小网格尺寸为0.05 m，总网格数约为900万. 计算域的入口边界设置为速度入口，出口为压力出口，洞体壁面与模型均采用无滑移壁面.

在模拟中，流体视为连续相，流场模拟选用Realizable k-ε湍流模型；依据模拟原型^[5-6]，入口风速为20 m/s，湍流强度与积分尺度分别为5%与0.01 m. 空气密度取1.396 kg/m³，黏度为1.61×10⁻⁵ ${\text{Pa}} \cdot {\text{s}}$. 稳态求解采用SIMPLE算法直至计算收敛，其中压力项与动量项均选用2阶迎风格式，时间项选用隐式1阶格式，湍流动能与湍流扩散率选用1阶迎风格式.

在稳态流场计算收敛后，在入口壁面建立颗粒工厂（2.1 m×1.0 m），直接顺来流方向引入模拟雪颗粒，进入DPM瞬态计算. 与原型一致，依据试验雪质量浓度平均值，换算颗粒工厂单位时间须入射的颗粒总量为100万. 模拟中，雪颗粒被视为球体^[4-14]，忽略了雪颗粒的自身破碎、形变. 颗粒粒径取均值250 μm，质量浓度为400 kg/m³，颗粒法向与切向恢复系数分别取0.01与0.2.

计算域的入口与出口均设置为逃逸边界，模型（楔体与固定架）设置为自定义捕捉边界，内嵌Trenker壁面覆雪判别准则^[4]. 雪颗粒引入时间与连续相时间步长一致（取10⁻⁴ s），最大时间步数设为5 000，每步的最大迭代次数为20，总计时长为1 s. 为了模拟湍流对雪颗粒运动轨迹的影响，计算中采用Fluent软件DPM颗粒随机游走模型.

如图5所示，楔体模型周边流场分布若干切片. 其中，切片1位于楔体模型A与B中间，受“窄道效应”的影响，两楔体模型间流速的增加显著，出现高流速区，平均流速约为25 m/s. 以楔体模型A为例，切片2穿过该模型的中心，楔体模型A迎风面的棱边处出现了明显的气流分离，分离后流速有所增加，对应的平均流速约为23.5 m/s. 楔体模型A背风面受模型阻塞的影响，在固定架下区域流速平均值仅为0.5 m/s.

在流场入口处引入雪颗粒，开展DPM瞬态计算. 如图6所示，入射的雪颗粒随流场运动，且颗粒场与流场的速度分布一致. 雪颗粒在两楔体模型间加速，在与模型壁面发生碰撞后减速. 楔体模型顶面、斜面法向的差异（差36°），造成雪颗粒在模型近表区域，受气流的分离与再附作用，发生明显的运动轨迹变化. 颗粒-壁面碰撞角 $\alpha $在楔体模型斜面不固定，呈现出两端 $\alpha $小、中部 $\alpha $大的特点. 相较而言，更多的雪颗粒在与楔体斜面的碰撞中，因不满足壁面覆雪准则^[4]，经反弹后再次返回流场.

以楔体模型A为例，提取各壁面单元的摩擦速度 $ {u _*} $，如图7(a)所示. 可知，该原型在20 m/s的试验风速条件下，模型顶面与斜面（均为迎风壁面） $ {u _*} $平均值分别约为0.85与0.70 m/s. 迎风壁面该参数的最大值达到1.05 m/s，位于靠近的两壁面的交界处. 由于模型顶面与固定架平面相连接（见图4），受壁面效应的影响，上述连接区域的边缘处摩擦速度较小.

通过阈值设定与像素提取可知，该楔体模型迎风壁面摩擦速度满足 $ {u _*} $<0.21 m/s的区域仅存在于顶面边缘内，如图7(b)所示. 在Trenker壁面覆雪判别准则^[4]中， ${u _{*{\rm{t}}}}$取值为0.21 m/s成立，楔体模型A颗粒捕捉（覆雪）的区域有且仅有极小的范围. 这与原型迎风全壁面均可以覆雪（见图3）的试验现象严重不符.

相较而言，蔡路等^[7-8]对参数 ${u _{*{\rm{t}}}}$的取值1.0 m/s更合理，可以近乎实现对模型迎风全壁面摩擦速度 ${u _*}$的包络，满足壁面均可以覆雪的条件.

从壁面覆雪判别准则^[4]可知，颗粒-壁面碰撞特性（ ${\alpha _{\rm{t}}}$）是决定壁面是否覆雪的又一要素，Trenker等^[4-9]对该参数的取值（45°）缺少合理验证. 在准则参数 ${u _{*{\rm{t}}}}$=1.0 m/s不变的前提下，对 ${\alpha _{\rm{t}}}$的取值开展梯度分析，讨论 ${\alpha _{\rm{t}}}$= 40~65°（每组间隔5°）情况下模型壁面的模拟覆雪差异.

以楔体模型A为例，如图8所示为在不同 ${\alpha _{\rm{t}}}$下模型迎风壁面的覆雪改变. 采用无量纲化的单元颗粒捕捉率 $\eta $来反映壁面覆雪特征，具体按下式计算：

当 ${\alpha _{\rm{t}}}$<50°时，楔体模型A迎风壁面的覆雪分布相似，主要集中在顶面与斜面尾端. 这是因为雪颗粒受模型近表气流的分离与再附作用，在斜面尾端与壁面发生碰撞， $\alpha $较小，易满足 $\alpha $< ${\alpha _{\rm{t}}}$的壁面覆雪准则要求. 随着 ${\alpha _{\rm{t}}}$的增大，楔体模型斜面的覆雪随着 $\alpha $分布变化，出现从两端向中间聚集的趋势. 直至 ${\alpha _{\rm{t}}}$大于60°，楔体模型A迎风壁面才出现较均匀的全壁面覆雪，与原型一致.

为了确定 ${\alpha _{\rm{t}}}$的取值，如图8(a)所示，提取楔体模型中轴线上的 $\eta $开展定量分析，如图9所示. 由于现行模拟方法不考虑壁面因覆雪累积而造成的外形改变，在确定的 ${\alpha _{\rm{t}}}$下，模拟所得的 $\eta $在模型壁面单元上保持稳定，可以反映出单元覆雪的累积规律. 当 ${\alpha _{\rm{t}}}$=45°时^[4]， $\eta $接近0，无覆雪累积. 当 ${\alpha _{\rm{t}}}$=60°时，楔体模型顶面与斜面的中部单元覆雪累积量比值约为1.51，十分接近原型试验（见图3）的覆雪结果.

在Trenker等^[4-8]的研究中，对壁面覆雪判别准则中 ${\alpha _{\rm{t}}}$的取值45° 较小，通过与原型的定性与定量对比，建议 ${\alpha _{\rm{t}}}$取60°.

为了平衡模拟计算代价与结果精确性的要求，基于DPM离散相模型配合自定义捕捉边界，对高速列车开展覆雪模拟研究，往往须对列车进行合理简化，主要包含转向架+棱台车体^[7-8]、转向架+截断车体^[4-5]2类.

选用我国现有运营的某型号高速列车，开展覆雪模拟分析. 建模时，对整车车体进行截断简化，仅对前车（含前、后2组转向架）建立模拟分析模型，如图10(a)所示.

模型计算域的尺寸为61W×31W×15W，高速列车模型距入口边界25W，W为标准车宽（约为3.1 m）. 空气相边界条件考虑车体与地面间的滑移，设置为移动（平动）壁面，轮对设置为移动（旋转）壁面. 列车运行速度取250 km/h，连续相空气参数的设置参考Liu等^[22]对哈大线沿线的实测值，质量浓度与黏度分别取1.453 kg/m³与1.57×10⁻⁵ ${\text{Pa}} \cdot {\text{s}}$. 雪颗粒物性参数按哈尔滨降雪测定值，质量浓度取320 kg/m³，粒径取均值500 μm. 其余连续相与离散相的边界设置及求解均与楔体模型的模拟一致，不再赘述.

模型采用非结构化网格，在前、后转向架区域进行网格加密、过渡处理，最小网格尺寸为1 mm. 近壁面采用标准壁面函数处理，壁面处y⁺为30~300，适应于Realizablek-ε模型下湍流核心区的充分发展. 为了验证模型网格的无关性，采用3种不同数量的非结构网格划分尺寸，分别记为粗网格、普通网格与细网格，对应体网格数为822万、1 058万与1 466万. 粗、细体网格划分设定值较普通网格放缩1.2倍，以普通网格划分为例，具体模型中轴面网格及壁面的加密如图10(b)所示.

考虑到模型附面的网格信息是直接影响 $ {u _*} $的重要因素， $ {u _*} $直接决定壁面覆雪的判别，因此选用转向架区域的 $ {u _*} $分布作为不同网格的检验参量. 以列车前转向架为例，结果如图11所示. 可以看出，在3种不同的网格尺寸下，检验区域 $ {u _*} $取值及其分布近乎相同，仅有粗网格在构架中心的分布有略微差别.

选用普通网格稳态流场的计算结果，开展DPM离散相模型分析. 颗粒工厂（引入雪颗粒）的设置与文献[9]一致，位于前转向架前侧. 单工厂尺寸为1.54 m×0.10 m，呈均匀点阵分布，垂直于列车车体平面，颗粒引入量设置为200万. 如图12所示为具体颗粒在流场中的运动迹线. 可以看出，受列车周边复杂流场的作用，气流将驱使雪颗粒从底部进入转向架区域. 流经端板分离后，将未实现壁面覆雪的颗粒从后部向两侧排出. 该现象（规律）与文献[14]类似的实测结果一致.

为了分析列车转向架区域流场的特征，截取整车截断模型中轴剖面顺风向速度分布，见图13(a). 可以看出，车体周边低流速区位于车体背风侧涡旋及车底近地表区域. 近地表低流速区会使该车底部位（如转向架）壁面的 $ {u _*} $减小，是导致该区域车体壁面覆雪的必要条件. 以流场发展更完全的后转向架为例，截取若干特征切片，对周边流场开展具体分析.

如图13(b)~(e)所示，高速列车运行时，由于列车转向架的构型复杂、构件繁多，内部空间流场发展表现出非定常、多涡旋的特点. 具体地，气流绕车体表面与地面空间加速流入转向架区域，该区域空间的增加及复杂的构件外形使得气流因分离发生方向的改变，形成扰流. 其中，扰流形成的低流速区域主要位于转向架半封闭“空腔”区域，如中心牵引装置内部（切片1）、构架自身或与其他构件（切片2~4）间隙区域，车底空间高流速区域集中在迎风侧构件（车轮、轮毂与齿轮减速箱）的立面. 这些区域的 $ {u _*} $相对较高，后轮毂处 $ {u _*} $最高可以达到约5.5 m/s（见图14）. 来流流经转向架后，受车体底部端板外形的影响，会在转向架后部形成明显的大尺度涡旋.

雪颗粒随流场从底部卷入上述低流速空间，作用在颗粒上的拖曳力改变，但颗粒运动来不及随流场变化而同步，在惯性作用下，与转向架各部件壁面发生碰撞. 经壁面覆雪准则（ $ {u _{*{\rm{t}}}} $=1.0 m/s与 $ {\alpha _{\rm{t}}} $ = 60°）的判别后，通过壁面单元颗粒捕捉率可以实现对列车转向架区域覆雪的评估，如图15所示.

结果表明，列车前、后转向架区域易出现覆雪的区域主要位于各构件背风侧壁面，尤以齿轮减速箱与构架横梁为甚，对应壁面的 $\eta $极值高达95%. 该分布规律与文献[8]的模拟结果高度一致. 由于列车前转向架更靠近车头鼻尖，车体前方来流相对受扰较小，二级悬挂上覆雪的区域以背风侧立面为主，平均壁面的 $\eta $约为78%，且下方的抗蛇行减振器无明显覆雪. 相比而言，列车后转向架周围来流因充分发展，致使该区域壁面覆雪分布更分散，对应的 $\eta $更小，易出现覆雪区域（构架横梁背风侧）壁面的 $\eta $极值仅为86%. 后转向架二级悬挂的内侧出现大范围覆雪，覆雪的面积约为前转向架该区域的3.5倍. 在后转向架抗蛇行减振器尾部出现小范围覆雪，对应壁面的平均 $\eta $为81%.

与文献[8]的分析一致，采用 $\eta $可以较好地评估转向架区域覆雪的发展，即由 $\eta $较高的转向架中部构架横梁处向前、后两侧发展，因此可以用于对高速列车减防雪评估、构架优化设计的工作. 对于模拟的真实转向架覆雪，须讨论现有壁面覆雪判别准则在DPM瞬态计算中的更新.

（1）建立壁面阈值摩擦速度与颗粒-壁面捕捉角的覆雪判别准则^[4]，考虑流场作用与颗粒-壁面碰撞两方面的特征，对雪颗粒的自定义捕捉，实现壁面覆雪模拟.

（2）通过控制变量参数分析法，对Paradot开展的楔体覆雪原型试验^[5-6]进行同边界模拟及准则参数讨论. 结果表明，针对原型雪颗粒及试验环境特征，对应模拟中壁面覆雪准则参数的取值偏小，建议 ${u _{*{\rm{t}}}}$与 ${\alpha _{\rm{t}}}$分别取为1.0 m/s与60°.

（3）利用改进的准则参数，定义无量纲化的单元颗粒捕捉率 $\eta $，对我国现有运营某型号高速列车运行时的转向架区域覆雪分布开展模拟评估. 结果表明，列车转向架易出现覆雪的区域主要位于各部件背风侧壁面及连接角落，尤以齿轮减速箱与构架横梁为甚，对应的 $\eta $极值高达95%. 相比而言，后转向架的周边流场因充分受扰，致使覆雪分布区域更分散， $\eta $相对更小.