高速列车在积雪轨道长时间运行时，被运行列车底部气流卷入的雪粒粘附和沉积到转向架表面上形成积雪，受转向架区域发热部件的影响，部分雪粒或积雪融化后形成的水滴运动到积雪面上会转化为积冰. 当转向架积雪结冰达到一定程度时，会阻碍悬挂系统的运动、挤压或损坏转向架上的零部件，减振系统、制动装置的功效会受到影响^[1].

转向架区域积雪与风致雪漂移运动有关，研究者们对地表大气边界层内风致雪漂移现象进行了大量的实验和数值仿真研究^[2-5]. Beyers等^[6-7]采用计算流体力学（computational fluid dynamics，CFD）方法和动网格技术，模拟地表方形建筑物周围的雪堆发展过程；周晅毅等^[8]利用欧拉多相流方法，分析航站楼屋盖表面的雪载荷分布；Sundsbo^[9]利用流体体积法（volume of fluid，VOF）模型分析气流通过台阶时的风雪流，研究导风板对积雪的影响；Thiis等^[10]利用CFD方法，研究三维柱状阵列防风林的积雪沉积效率；Tominaga等^[11]对立方体建筑模型周围的雪堆进行CFD预测. 计算流体动力学方法已经广泛应用于模拟建筑物周围的风雪流^[12]. 上述研究中，对雪相的处理分为欧拉方法和拉格朗日方法两类；相对于欧拉方法，拉格朗日方法可以更加直观地从单个颗粒的角度对雪粒的运动机理进行仿真研究^[13].

与大气边界层内风致雪漂移运动不同，高速列车车底雪粒漂移主要由列车风引起. Kwon等^[14]的研究表明，在距离轨道地面0.18 m的高度，列车引起的列车风速度可以达到车速的40%；Shishido等^[15]在降雪风洞中研究列车裙板侧导流板对转向架区域积雪的影响；Ido等^[16]在低温风洞中研究列车裙板侧三角菱形扰流板对转向架区域积雪的影响；Xie等^[17]利用RNG $k - \varepsilon $湍流模型和离散相模型（discrete phase model，DPM），研究高速列车转向架区域的积雪特性；Wang等^[18-19]采用DPM，研究转向架舱端墙斜度、雪粒密度、雪粒直径对转向架积雪量的影响；Paradot等^[20]利用与来流成不同角度的平板的风洞积雪试验数据，建立基于雪粒入射角度的雪粒捕获准则，并用于高速列车转向架积雪模拟；Gao等^[21-22]利用CFD方法，分析转向架的积雪分布，用木屑模拟雪粒进行转向架积雪试验；Trenker等^[23]采用DPM模型和自定义雪粒捕获边界，分析列车底部设备上的积雪速率. 目前的列车风雪流研究侧重于转向架区域雪粒运动轨迹的分析，高速列车底部气流引起的雪粒输运与堆积现象的研究较少.

由于缺乏可以通用的雪粒堆积模型，目前转向架区域的积雪量的计算多采用壁面全捕捉^{[17-19, 21-22]}，这与实际情况存在一定差异. 本文建立基于临界捕获角度和临界剪切风速的雪粒捕获模型，用于处理雪粒与壁面间的碰撞行为及计算壁面积雪量，采用DPM方法对国内某型动车组拖车转向架区域的风雪两相流场进行数值模拟，分析转向架区域的雪粒运动和沉积特性.

由于雪粒的运动模拟不适用相似性准则^[17]，不宜采用缩比模型模拟风-雪两相流，若采用3节编组的全尺寸列车模型，DPM方法需要追踪超千万级的粒子轨迹，会消耗非常多的计算资源. 考虑到远离车头、车尾的转向架区域流动受上游流场及车体顶部外形的影响较小^[21-22]，为了节省计算资源，仅取转向架局部区域进行计算. 以某型动车组中间车二位端拖车转向架为研究对象，将车体部分简化为车体罩（长L=10 m），侧面外形与整车相同，前后端面与水平面成30°倾斜角，几何模型如图1所示. 转向架外流场计算区域的长度为7L、高度为3L、宽度为5L，转向架罩前端距离入口2L、底部距离轨面0.2 m. 计算域前端和后端分别设置为速度入口边界和压力出口，顶部和两侧面为对称边界，地面和轨道设置为移动壁面.

气相采用RANS（Reynolds averaged Navier Stokes）方程，湍流模型采用Realizable $k - \varepsilon $模型^[24]. 该湍流模型对于旋转流动、强逆压梯度的边界层流动、流动分离和二次流有很好的表现，风洞测试结果^[21]表明，该湍流模型用于计算转向架区域流场的准确度较高. 雪粒运动轨迹采用DPM求解，该模型中流体的运动用连续介质的N-S方程进行描述. 颗粒的运动通过在拉格朗日框架下分析每个粒子的运动轨道进行描述，粒子的速度及位移的更新由邻近颗粒对它的碰撞过程及周围流体对它的悬浮过程来确定，颗粒对流体的瞬时作用反映在N-S（Navier-Stokes）方程求解过程中不断修正的源项上. 假设雪粒为球体，不考虑雪粒间碰撞的影响. 雪粒的力平衡方程如下：

(1) $\frac{{{\rm d}{{ u}_{\rm p}}}}{{{\rm d}t}} = {F_{\rm D}}\left( { u - {{ u}_{\rm p}}} \right) + \frac{{ g\left( {{\rho _{\rm p}} - \rho } \right)}}{{{\rho _{\rm p}}}} + F.$

式中： ${F_{\rm D}}\left( { u - {{ u}_{\rm p}}} \right)$为粒子受到的拖曳力； ${F_{\rm D}}$为粒子松弛时间的倒数，定义为

(2) ${F_{\rm D}} = \frac{{18\mu }}{{{\rho _{\rm p}}d_{\rm p}^2}}\frac{{{C_{\rm D}}{{Re}} }}{{24}},$

(3) ${{Re}} = \frac{{\rho {d_{\rm p}}\left| {{{ u}_{\rm p}} - u} \right|}}{\mu }.$

粒子采用光滑球模型^[25]，阻力系数 ${C_{\rm D}}$定义为

(4) ${C_{\rm D}} = {a_1} + \frac{{{a_2}}}{{{{Re}} }} + \frac{{{a_3}}}{{{{{{Re}} }^2}}}.$

式中： $ u$为空气速度； ${ u_{\rm p}}$为雪粒速度； $\mu $为空气流体动力黏度； $\rho $为空气密度； ${\rho _{\rm p}}$为雪粒密度； ${d_{\rm p}}$为雪粒直径； ${a_1}$、 ${a_2}$、 ${a_3}$为与粒子雷诺数Re范围有关的常数，取值见文献[25].

研究雪粒在转向架区域的悬移运动以及雪粒与转向架壁面的撞击行为，不考虑车底风雪流的起动过程. 由于数值模拟采用风洞方式，以列车为参照系，实际运行时离开积雪路面的雪粒速度相对于列车的运动速度为小量，雪粒相对于列车的运动速度接近车速，故数值计算时用固定在列车底部的垂向喷射面向列车底部区域水平喷射初始速度为车速（入口风速）的粒子. 雪粒喷射面布置在转向架罩下部，位置如图2所示，高度为0.4 m，宽度为4.0 m. 粒子喷射器参数如下：初始速度为（V_t，0，0），质量流率为0.1 g/s，雪粒直径取200 μm，密度取150 kg/m³.

雪粒是否粘附在壁面上取决于雪粒与壁面以及气流的相互作用，表现为雪粒与壁面发生碰撞时的接触力和附近气流对其产生的流体动力. 雪粒在壁面上是否粘附或反弹，这一行为取决于雪粒的碰撞能力（反弹能力）和壁面的吸附能力；雪粒在壁面上能否稳定沉积，这一行为取决于雪粒的临界剪切速度和气流的摩擦风速. 天然雪粒形状不规则、结构松散，且不同类型的雪粒的力学性质受雪粒形状、雪粒黏性、雪晶体间的黏结特性等因素的影响较大，很难用统一的反弹恢复系数或击溅函数描述雪粒与壁面的碰撞行为. 参考文献[23]的雪粒捕获准则，建立基于临界捕获角度和临界摩擦风速的雪粒沉积模型：当 ${\alpha _{\rm p}}$< ${\alpha _{\rm{trap}}}$时，雪粒粘附于表面，否则发生反弹；若雪粒撞击壁面处的摩擦风速 ${u_ * }$大于临界摩擦风速 ${u_{ * {\rm t}}}$，则雪粒不会发生稳定沉积，粒子会被剥离表面；雪粒与壁面的碰撞结果如图3所示.

Beyers等^[26]在南极地区测得的陈年的风硬化积雪场的临界摩擦风速为0.28 m/s，指出新的积雪场的取值为0.20~0.28 m/s. Clifton等^[27]的研究表明，雪密度和尺寸对 ${u_{ *{\rm t}}}$的影响大于环境温度和湿度的影响，测得的 ${u_{ *{\rm t}}}$为0.27~0.69 m/s. Kind^[28]的研究表明，新降的干燥雪粒的 ${u_{ *{\rm t}}}$为0.07~0.25 m/s，经过一段时间风侵蚀后，临界摩擦风速为0.25~1.0 m/s. 考虑到高速列车运行速度高，雪粒撞击壁面速度较大，经过长时间的风雪侵蚀后，积雪表面临界摩擦风速会较大. 本文 ${u_{ *{\rm t}}}$取Kind测得的临界摩擦风速中的极大值1.0 m/s. 计算时，临界捕获角度取π/4 rad^[23]，临界摩擦风速取1.0 m/s，假定被气流剪切走的雪粒获得的速度为临界摩擦风速.

为了分析转向架壁面上的雪粒沉积特征，在Fluent软件中定义4个自定义存储器（user defined memory，UDM），存储壁面上的雪粒撞击数、雪粒入射质量、雪粒堆积量、雪粒黏附率. 定义单位面积上撞击壁面的粒子质量为雪粒入射质量，定义单位面积上粘附于表面的粒子质量为雪粒堆积量，定义粘附于表面的粒子质量与撞击表面的粒子质量之比为黏附率. 基于Fluent软件计算气相流场和粒子轨迹，雪粒与壁面的碰撞行为采用自定义边界条件（user defined functions，UDF）实现，壁面上相关统计量的计算结果存入UDM中. 雪粒与壁面碰撞行为处理流程如图4所示. 1）计算 ${u_ * }$、 ${\alpha _{\rm p}}$；2）判断 ${\alpha _{\rm p}}$是否小于 ${\alpha _{\rm{trap}}}$：若是，则捕获粒子等待进一步判断，否则计算新的反弹速度，并激活雪粒轨迹；3）判断被捕获雪粒处的 ${u_ * }$是否小于 ${u_{ *{\rm t}}}$：若是，则雪粒判定为堆积、计算壁面处堆积量；否则粒子被气流剪切走，假定雪粒获得与摩擦风速同等的风速，并激活雪粒轨迹；4）进入下一时间步，连续相和粒子相时间步迭代完成后循环执行1）、2）、3）.

流体区域离散采用四面体网格，对简化车体周围、尾流区域进行2次网格加密. 近壁区采用壁面函数处理，壁面附近划分10层边界层网格，每层的厚度增长率为1.2，使壁面处的y⁺为30~300. 设置3套（粗网格、中等网格、细网格）不同密度的网格进行独立性检验，体网格数量分别为1.8×10⁷、2.4×10⁷、3.4×10⁷. 粗网格和细网格的尺寸相对于基准网格分别扩大和缩小了1.2倍，但边界层尺寸设置相同. 取转向架纵向对称面上距轨面0.1 m高的水平中心线上的压力系数C_p作为不同网格的检验量，结果如图5所示. 通过对比取中等密度网格进行后续的计算，网格如图6所示.

连续相采用Realizable $k - \varepsilon $湍流模型、相关参数采用软件推荐值. 雪粒采用DPM，雪粒-壁面边界参数如下：捕获角度为45°，临界摩擦风速为1.0 m/s，法向和切向反弹系数分别为0.01、0.2^{[19, 22]}. 计算模型边界条件设置如表1所示，其中转向架除轮对外粒子相的边界条件采用自定义边界条件（UDF）. 列车运行速度V_t=300 km/h，空气参数参考哈大线极端低温环境（−30 °C）下的实测值^[29]： $\mu $=1.57×10⁻⁵ Pa·s， $\rho $=1.453 kg/m³. 先计算稳态流场，收敛后再开始瞬态计算，时间步长为10⁻⁴ s，每个时间步内迭代20次，总求解时间为1.0 s. 采用SIMPLE（semi-implicit method for pressure linked equations）算法求解，压力项采用Standard格式，动量项使用2阶迎风格式，湍流动能和湍流扩散率采用1阶迎风格式，时间项采用隐式1阶格式.

粒子采用非稳态粒子追踪方法，雪粒子注入时间与连续相时间步同步（ $\Delta t = 10^{ - 4}\;{\rm s}$），选择随机游走模型（random walk model，RWM）模拟湍流脉动对雪粒运动轨迹的影响. 粒子轨迹最大计算时间步数（maximum number of steps）设为5 000，满足粒子流经计算域所需的总时间步数. 粒子运动方程积分时间步控制参数 $\lambda $（step length factor）设为10，使其满足 $\lambda {{N}} \approx $N_max，其中N为雪粒沿流动方向经过的网格数，N_max为最大时间步.

为了分析转向架区域气流方向，在转向架区域作4个纵向速度切片，切片位置如图7所示，切片上流线如图8所示. 图中，u_x为x方向上的速度. 由图8可以看出，受旋转轮对影响，车体底部气流主要由转向架前轮对和后轮对后方向上扬进入转向架内部区域，并在两轮对周围形成回旋气流；转向架后下部为主要的气流进入区域、前部为主要的气流流出区域；进入转向架的气流速度较高，出转向架区域的气流速度较低，转向架区域气流速度分布整体上呈现出后下部速度较高、前上部速度较低的特征. 低速气流本身挟带雪粒运动的能力较弱，高速气流较强，由于雪粒良好的随流特性，转向架区域气流的雪粒输运能力呈现出进强出弱的特征.

雪粒进入转向架区域的过程如图9所示. 可以看出：由车底气流挟带的雪粒在转向架区域底部区域减速，并由两轮对后部向上翻转进入转向架区域，部分雪粒被轮对、刹车盘边缘带动绕车轴旋转；车底下层的距离转向架部件较远的雪粒会被车底气流直接挟带离开转向架区域. t=1.0 s时转向架区域的雪粒分布如图10所示，雪粒子采用粒子停留时间t_p（particle residence time）染色. 可以看出，转向架区域大部分雪粒停留时间小于0.6 s，前端墙处的雪粒停留时间较长，表明这部分雪粒在转向架区域运动了较长时间才到达前端墙区域.

临界摩擦风速是风可挟带雪粒运动的阈值^[4]，若壁面处湍流流动摩擦风速大于雪粒临界摩擦风速，积雪表面会发生侵蚀，则沉积的雪粒将会被气流吹走. Kind的研究表明，经过一段时间风侵蚀后的临界摩擦风速为0.25~1.00 m/s^[28]. 当列车速度为300 km/h时，转向架表面的摩擦风速如图11所示. 可以看出，转向架大部分区域表面处的摩擦风速小于临界摩擦风速1.0 m/s，摩擦风速小的表面雪粒易发生堆积；旋转部件轮对、刹车盘及转向架底部的部分气流冲击区域摩擦风速很大，不易发生积雪.

0.2 s时间内转向架表面上的雪粒子入射质量M_i分布如图12所示（云图显示范围为0~0.2 g/m²，全局范围为0~735.2 g/m²）. 由图12可以看出：转向架底部受来流直接冲击的部件表面雪粒入射质量很大，其中构架横梁迎风面雪粒的撞击量最大，0.2 s内入射质量超过0.2 g/m²；转向架后轮对（轮对2）区域受向上翻转气流冲击的部件表面雪粒入射质量很大，其中摇枕后壁面、牵引拉杆后部面、后轮对中间制动夹钳迎风面雪粒撞击量最大，0.2 s内入射质量超过0.2 g/m²；侧面的抗蛇行减振器、后部轮对的轴箱悬挂等部件壁面上雪粒入射质量较大，为0.1~0.2 g/m². 构架底部、抗蛇行减振器、后轮对中间制动夹钳、牵引拉杆、抗侧滚扭杆区域为雪粒子集中撞击部位.

0.2 s时间内不同部件接收到的撞击雪粒总数如图13所示. 图中，N_i为总雪粒撞击数. 可知，轮对、构架上的撞击雪粒数超过10⁶；制动夹钳、刹车盘、前轮对轴箱悬挂、抗蛇行减振器上的撞击雪粒数超过10⁵；摇枕、后轮对轴箱悬挂上的撞击雪粒数超过10⁴；踏面清扫装置、空簧、抗侧滚扭杆上的撞击雪粒数超过10³. 构架底部、制动夹钳下部等迎风区域为雪粒子集中撞击部位，底部区域粒子撞击数比上部区域高1、2个数量级.

转向架表面上0.2 s内雪粒堆积质量M_a分布如图14所示（云图显示范围为0~0.1 g/m²，全局范围为0~22.9 g/m²）. 可以看出，转向架底部受来流直接冲击的部件表面雪粒堆积质量很大，其中构架横梁迎风面雪粒堆积量最大，0.2 s内堆积质量超过0.1 g/m²；后轮对（轮对2）区域受向上翻转气流冲击的部件表面雪粒堆积质量很大，其中摇枕后壁面、牵引拉杆后部面、后轮对中间制动夹钳迎风面雪粒堆积量最大，0.2 s内堆积质量超过0.1 g/m²；侧面的抗蛇行减振器、后部轮对的轴箱悬挂等部件壁面上雪粒堆积量较大，为0~0.05 g/m². 构架底部、抗蛇行减振器、后轮对中间制动夹钳、牵引拉杆、抗侧滚扭杆区域为雪粒易堆积部位，这些部件为关键悬挂及制动部件.

转向架各部件上的总积雪量与逃逸量（总入射量−总堆积量）的占比如图15所示. 图中， $\overline M_{\rm a}$为平均堆积量， γ_a、γ_e分别为积雪量和逃逸量的占比. 可以看出，撞击在一位端轴箱悬挂组件上的雪粒约50%发生了沉积；空气弹簧上有30%的撞击雪粒发生了沉积；抗侧滚扭杆组件上有20%的撞击雪粒发生了沉积；踏面清扫装置上有13%的撞击雪粒发生了沉积；摇枕上有12%的撞击雪粒发生了沉积；抗蛇行减振器上有12%的撞击雪粒发生了沉积；其他部件的雪粒沉积比例均小于10%. 各部件上的平均堆积量（总堆积量/总面积）按大小排序依次为：牵引拉杆、构架、摇枕、制动夹钳2、抗侧滚扭杆、制动夹钳1、轴箱悬挂2、横向减振器、轴箱悬挂1、空气弹簧、抗蛇行减振器、踏面清扫装置. 其中，牵引拉杆、构架、摇枕、制动夹钳2、抗侧滚扭杆、制动夹钳1的平均堆积量大于0.01 g/m²（0.2 s时间内）.

转向架表面雪粒黏附率R_a分布如图16所示. 可以看出：前轮对区域构架顶部上表面雪粒撞击区域黏附率接近1.0，后部轮对区域的构架上表面黏附率小于0.5，造成差别的主要原因是列车底部气流主要由转向架后部向上翻卷进入转向架区域，转向架前轮对所在区域的气流速度比后轮对区域小，雪粒沉积在前轮对区域顶部壁面的概率大. 构架横梁与弓形侧梁下部弯角处、抗蛇行减振器安装座角落处、摇枕拐角处的最大黏附率能够达到1.0左右. 牵引拉杆、后部制动组件（夹钳2）、抗侧滚扭杆、空气弹簧及构架横梁等受上扬翻转气流冲击的部件表面，最大黏附率能够达到1.0左右. 黏附率最小的是转向架部件底部面，因其表面与气流夹角较小，雪粒不易堆积.

总体而言，转向架后部区域垂向面黏附率高，前部区域水平面黏附率高，角落区域黏附率高. 无论是积雪量，还是黏附率，都是以转向架中部横梁区域为最大. 该区域集中安装了制动组件、牵引拉杆以及除一系悬挂外的所有悬挂部件，若积雪转化为积冰，则会严重阻碍悬挂部件运动，威胁行车安全. 从仿真结果来看，转向架构架横梁区域由于积雪量大且雪粒黏附率高，积雪将会以转向架中部区域开始发展，且后侧积雪速度快于前侧. 实际中转向架的积雪情况如图17所示. 可以看出，转向架中部区域以及制动夹钳周围积雪严重，这与前面的分析结果一致.

不同捕获角度对应的各部件的平均入射质量（总入射质量/部件总表面积）、平均堆积量（总堆积量/部件总表面积）、平均黏附率（总堆积量/总入射质量）如图18所示. 图中， $\overline M_{\rm i} $为平均入射量， $\overline R_{\rm a} $为平均黏附率. 可以看出，各部件平均积雪量由高到低依次为：牵引拉杆、构架、摇枕、夹钳2、抗侧滚扭杆、夹钳1、横向减振器、轴箱2、轴箱1、空气弹簧、抗蛇行减振器、踏面清扫，其中牵引拉杆、构架、摇枕、夹钳2的平均积雪量比其他部件高出约1倍. 转向架后部比前部积雪量大，其中后制动夹钳的平均积雪量比前制动夹钳高出约1倍. 当捕获角度增大时，各部件的雪粒平均入射量变化趋势不同：其中牵引拉杆、构架、摇枕、抗蛇行减振器、空气弹簧等部件的雪粒总入射量随着捕获角度的增大而减小；夹钳1、轴箱2、横向减振器、抗侧滚扭杆等部件的雪粒总入射量随着捕获角度的增大而增大；夹钳2、轴箱1、踏面清扫等部件的雪粒总入射量随着捕获角度的变化不明显；当捕获角度增大时，各部件的积雪率基本呈增大趋势，其中轴箱1和空气弹簧的积雪率对捕获角度的变化最敏感，增幅分别为24.8%和15.9%；当捕获角度增大时，各部件的平均堆积量基本呈增大趋势，但增加幅度不大. 当捕获角度为30°~60°时，捕获角度的变化对部件总的积雪量计算结果影响较小. 不同捕获角度下的转向架表面的雪粒堆积量分布如图19所示. 可以看出，捕获角度的变化对积雪区域的分布影响较小.

（1）转向架区域气流整体上呈现出后下部速度较高、前上部速度较低的特征；转向架后下部为主要的气流进入区域、前部为主要的气流流出区域；进入转向架的气流速度较高，流出转向架区域的气流速度较低，转向架区域气流在速率上呈现出进强出弱的特征.

（2）转向架的构架底部、抗蛇行减振器、制动夹钳、牵引拉杆、抗侧滚扭杆区域为雪粒易堆积部位. 转向架后部区域垂向面黏附率高，前部区域水平面黏附率高，角落区域黏附率高. 无论是积雪量，还是黏附率，都是以转向架中部横梁区域为最大，积雪将会以转向架中部区域开始发展，且后侧积雪速度快于前侧.

（3）转向架各部件平均堆积量由高到低依次为：牵引拉杆、构架、摇枕、二轴制动夹钳、抗侧滚扭杆、一轴制动夹钳、横向减振器、二轴轴箱、一轴轴箱、空气弹簧、抗蛇行减振器、踏面清扫装置，其中牵引拉杆、构架、摇枕、二轴制动夹钳的平均积雪量比其他部件高出约1倍，二轴制动夹钳的平均积雪量比一轴制动夹钳高出约1倍.

（4）当临界捕获角度为30~60°时，捕获角度的变化对部件总的积雪量的计算结果影响较小. 不同类型的雪粒力学性质受到雪粒形状、雪粒黏性、雪晶体间的黏结特性等因素的影响，差异较大. 本文在处理入射角度大于临界捕获角度的雪粒时，采用统一的反弹恢复系数处理雪粒与壁面的碰撞行为，这会与实际存在一定差异，如何更加真实地描述雪粒与壁面的碰撞行为需要进一步的研究.