铁路缓和式路堤的风沙流场数值模拟研究

徐建刚，黄 宁，石广田，张小安

(1. 兰州交通大学 铁道技术学院，兰州 730070；2. 兰州交通大学 甘肃省轨道交通服役环境与智能运维重点实验室，兰州 730070；3. 兰州大学西部灾害与环境力学教育部重点实验室，兰州 730000；4. 兰州交通大学 机电工程学院，兰州 730070)

随着我国“一带一路”经济带建设发展和“八纵八横”铁路网逐步实现，西部地区的铁路线路不断增多，如兰新铁路、青藏铁路及太中银铁路等，这些穿越戈壁、荒漠及沙漠地区的铁路经常遭受风沙危害，对列车的运行安全造成影响[1-4].

铁路风沙灾害的形成原因主要包括沙物质来源、风沙动力环境及施工设计不合理等[5-8]，因此铁路沿线设置了多种防风治沙措施保障列车安全运行[9].挡沙墙和挡风墙是应用最广泛的两种机械防风沙措施，学者们采用理论分析和风洞实验等方法，对铁路沿线挡沙墙、HDPE防沙网、挡风墙的防风沙机理和效果进行了大量研究[10-18].

铁路防风沙措施的最根本目的是减少路堤上轨道积沙，降低沙埋轨道等危害.因此，学者们采用数值模拟、风洞试验及野外实测等手段，对铁路路堤及各类沙障周围风沙流特性进行了广泛研究.石龙等[19]研究了风沙两相流对铁路路堤的响应规律，结果表明路堤顶面积沙量呈中间多两边少的“正态形”分布，且风速增大路堤迎风坡积沙量减小，背风坡积沙无明显变化；张军平等[20]探讨了兰新铁路戈壁地区路基周围风沙流运动特征，并给出了防风挡沙墙的合理高度应满足的若干条件；李晓军[21]研究了风沙两相流对铁路路堤及防风挡沙墙的响应规律，对比分析了不同类型挡风墙背风侧流场特征和线路积沙形态，并优化了挡风墙设计参数；邓腾飞[22]通过建立详实的有砟、无砟轨道模型，研究了有砟、无砟轨道周围的风沙流运动规律和沙粒沉积过程，得出沙粒在无砟轨道通过性更优的结论；秦旗[23]依据阿联酋地区的实际风沙条件，采用CFD技术模拟并给出了该线路路基周围多种风沙防治措施设计参数的适宜范围.刘刚等[24]提出在沙漠公路风沙病害治理中可采用流线型、圆弧化截面路基作为一种疏导风沙的措施，减少风沙运动阻力，促使风沙快速通过路堤.

综上所述，关于铁路风沙灾害防治的研究目前主要集中在各种防风沙措施对路堤周围流场、沙粒沉积及其防护效益的影响等方面.然而，铁路路堤自身的结构形式对风沙防护效果的影响未引起人们的足够关注，文献[24]中虽然提出流线型、圆弧化路基断面可作为沙漠公路的输沙措施，但目前仍没有相关验证，而在铁路沙害的治理措施中也很少有学者提及.因此，本文基于双欧拉流体模型，对比分析了普通路堤与缓和式路堤周围的流场分布特征和沙粒沉积状态，得出有利于减少沙粒沉积的路堤结构，以期为今后铁路防沙工程提供参考.

1 理论方程

本文主要采用欧拉双流体模型研究近地表风沙运动状态[25-26]，该模型视计算域内的流固两相物质为互相贯穿的连续介质，两相体积分数之和为 1 ，气相和沙相分别由各自的质量和动量守恒方程控制，且存在相互作用力.

αg+αs=1.

(1)

式中，αg、αs分别为气相、沙相的体积分数.

气相质量守恒方程：

(2)

沙相质量守恒方程：

(3)

气相动量守恒方程：

(4)

沙相动量守恒方程：

(5)

气固相间作用力：

(6)

式中，ρg、ρs分别为气、沙相的密度；t为时间；νg、νs分别为气、沙相的速度矢量；Ug、Us分别为气、沙相速度；τg、τs分别为气、沙相的表面应力张量；ps为沙相固体压力；g为重力加速度；fsg为气固相间作用力；CD为阻力系数；Ur为气固两相间的相对速度；d为沙粒的当量直径.

流体在近地面或者地形不规则处会出现运动复杂的涡流，多个涡流叠加起来就会形成湍流.气固两相流满足标准k-ε模型[27]的湍动输运方程：

(7)

(8)

式中，kg、μt、εg分别为湍动能、湍动粘度、湍动耗散率；Ugi、Ugj为速度在x、y方向上的分量；σk、σs为湍动能、湍动耗散率对应的普朗特数；Gkg为平均速度梯度引起的湍动能产生；Gb为浮力影响引起的湍动能产生；Ym为可压缩湍流的耗散率比；C1ε、C2ε为常数.

2 数值模拟

2.1 路堤模型

普通铁路路堤的两侧边坡为斜线坡，而缓和式路堤的边坡还可设计为下凹坡或上凸坡.其中，左侧下凹坡、上凸坡是以普通路堤的P1为起点，P2为端点，包含角分别为30°、-30°的圆弧；右侧下凹坡是以P3为起点，P4为端点，包含角分别为30°的圆弧.采用a、b、c分别代表斜线坡、下凹坡及上凸坡，以两两组合的方式能够构成9种路堤结构.其中，a-a路堤是普通路堤，再选取三组有代表性的缓和式路堤，分别为左侧斜线坡右侧下凹破的a-b路堤、两侧下凹坡的b-b路堤及左侧上凸坡右侧下凹破的c-b路堤，图1给出了四种路堤结构示意图.采用CAD绘制a-a路堤的几何模型及计算域示意图(图2).路堤的高度、顶部及底部宽度分别为5 m、8 m、25.5 m，计算域尺寸为50 m×155.5 m，计算域左侧入口到路堤迎风坡坡脚的距离为50 m，保证风沙流在此范围内充分发展.其他三组缓和式路堤仅两侧边坡形式不同，路堤及计算域的基本尺寸保持不变.

图1 四种路堤示意图Fig.1 Schematic diagram of four kinds of embankments

图2 a-a路堤的几何模型及计算域示意图Fig.2 Schematic diagram of geometric model of a-a embankment and simulation domain

利用ICEM软件的Patch Dependent方法对a-a路堤几何模型及计算域进行网格划分，划分成四边形为主(Quda Dominant)的结构化网格，并对路堤顶部区域进行局部加密，网格数约为150 000个，a-a路堤的网格划分结果如图3所示.其它三种路堤采用同样的方法进行网格划分.

图3 a-a路堤的网格划分结果Fig.3 Meshing result of a-a embankment

2.2 边界条件及求解方法

利用FLUENT15.0软件对风沙流进行模拟，分别导入四种路堤有限元模型，设置流场左侧边界为速度入口(VELOCITY_INLET)，右侧边界为自由流出口(OUT_FLOW)，上边界采用对称边界条件(SYMMETRY)，地面、路堤表面设置为固体无滑移边界(WALL)，下壁面粗糙高度设置为0.000 6 m[28].利用用户自定义函数(UDF)设置流场入口气流边界条件为典型风速廓线流：μ=μ*/κln(z/z0)，其中μ为高度z处的水平风速，μ*为摩阻风速，κ为冯·卡门系数，z为高度，z0为粗糙长度.流场入口处的湍流强度设置为5%，入口处沙相体积分数设置为0.01.

基于有限体积法进行离散，采用QUICK离散格式，利用SIMPLE压力-速度耦合修正算法，选择基于压力法的求解器，时间步长设置为0.01 s，收敛标准为残差小于1.0×10-6.风沙流场模拟中其它关键参数设置如表1所示.

表1 其它关键参数Tab.1 Other key parameters

2.3 流场验证

为了验证数值模拟中流场设置的可靠性，本文建立与文献[15]中风洞实验段纵断面尺寸相同的流场计算域：22 m(长)×1.45 m(高)，计算摩阻风速μ*为0.639 5m/s时，流场中近地表的风速廓线，并与文献[15]中风洞试验所得风速廓线结果进行对比，如图4所示.由图4可知，本文计算的风速廓线与风洞试验结果一致性较好，验证了数值模拟中流场设置合理、可靠，后续将采用相同的流场设置对各种路堤进行模拟计算.

图4 数值模拟与风洞试验的风速廓线对比Fig.4 Comparison of velocity profile between numerical simulation and wind tunnel test

3 结果分析

3.1 路堤周围流场特征

在风的剪切作用下地面沙粒会产生滚动、跳跃甚至飞入空中的现象，风沙动力学上对应称之为沙粒的蠕移、跃移与悬移运动[27].风是沙粒运动的动力源，风速对沙粒的运动起着决定性的作用.因此，本文利用Fluent模拟了摩阻风速为0.892 m/s，计算时间为20 s时，四种路堤周围的流场分布，并通过Tecplot软件处理后得到四种路堤的速度云图(图5).由图5(a)可知，路堤对于由远及近的风沙流来说是一个“减速带”，明显减弱了近地面、路堤迎风侧的风速，导致形成减速区；路堤迎风侧肩部对路堤顶部流场有明显加速作用形成加速区，同时路堤顶部较低区域由于涡流的出现形成低速区；路堤的遮蔽作用使背风侧较大区域形成回流区.图 5(b)～(d)中a-b、b-b及c-b三种路堤的流场分布与a-a模型类似，不再赘述.

图5 四种路堤周围的速度云图(单位：m/s)Fig.5 Velocity contours of four embankments (Unit:m/s)

铁路沿线风沙灾害中最直接影响列车正常运行是轨道沙埋，主要发生在路堤顶部轨道附近.因此，为进一步探明四种路堤在顶部轨道周围的流场分布，本文提取FLUENT计算结果，采用Origin软件绘制距路堤顶部0.1 m、0.5 m高度处水平风速沿程变化曲线，如图6所示.图中A、B、C、D分别表示风速在水平方向到达路堤迎风坡坡脚、路肩、顶部轨道及背风侧回流区等位置.图6(a)中四种路堤的风速变化过程如下：临近A点之前风速逐渐降低，随后气流受到路堤迎风坡的阻碍，风速在B点之前回升；气流在C点附近出现涡流，风速急剧下降至最小值，CD两点之间的风速逐渐上升，最后风速基本达到稳定.从整体可知，四种路堤的风速按照c-b、a-b、a-a及b-b的顺序依次减小，说明b-b路堤对气流的减弱效果最优.

此外，由图6(a)四种路堤在顶部轨道处的最小风速值可知，a-a路堤最小风速值为-4.0 m/s，而b-b路堤是a-a路堤的2.3倍，值为-9.1 m/s，气流在此处形成强涡流，且风速已达到沙粒的起动风速.因此，b-b路堤有助于气流在顶部轨道形成强涡流，使沙粒再次起动减少轨道积沙.图6(b)中距各种路堤顶部0.5 m处风速的变化规律与图6(a)基本类似，由于距路堤顶部高度的增加，风速值均大于零.

图6 四种路堤的水平风速变化Fig.6 Horizontal velocity variation of four embankments

图7给出了b-b路堤不同高度的水平风速沿程变化曲线，显示路堤背风侧有明显的回流区.图中E、F分别表示风速在水平方向达到路堤背风坡中部和回流区尾部位置.由图7(a)可知，距地面高度为4.0 m时，风速较大，回流区较小；而距地面高度≤2.0 m时，E、F之间区域形成较大回流区，回流区长度约45 m，风速降到最低，随后逐渐增加.在距地面0.5 m时，回流区中最小值风速为-11 m/s，有利于沙粒再次起动向远场输送.由图7(b)可知，距离路堤顶面高度≤0.3 m时，路堤顶部的风速降至最小出现负值；而距离路堤顶面高度≥0.5 m时，路堤顶部风速均为正值，风速变化规律与图6(a)一致.

图7 b-b路堤不同高度处的水平风速变化Fig.7 Horizontal velocity variation of b-b embankment at different heights

3.2 路堤周围积沙状态

路堤影响风沙流运动过程，路堤周围风速降低，导致沙粒沉积.积沙分布能够直观地反映路堤周围的积沙状态.本文给出了摩阻风速为0.892 m/s，计算时间为20 s时，四种路堤的积沙分布，通过Tecplot软件绘制如图8所示.

图8中路堤周围均出现沙粒沉积，主要分布在路堤迎风坡、背风坡坡脚及顶部轨道处.图8(a)中a-a路堤顶部沙粒主要沉积在2轨道右内侧，沙粒沉积严重，甚至达到钢轨的高度，将会对列车的安全运行产生较大影响；图8(b)中a-b路堤顶部沙粒在1、2轨道左外侧均有较少沉积，对列车行驶影响较小；图8(c)中b-b路堤顶部几乎没有沙粒沉积，不影响列车行驶；c-b路堤与a-b路堤顶部的积沙状态基本一致.由此可见，b-b路堤顶部积沙分布最少.

图8 不同路堤周围的积沙分布Fig.8 Sand deposition distribution around different embankments

此外，分析路堤两侧边坡形式不同对顶部轨道积沙的影响.a-a、a-b路堤的左边坡均设置为斜线坡，右边坡为下凹坡会使2轨道处的积沙向1轨道转移；a-b、b-b路堤的右边坡均设置为下凹坡，左边坡为下凹坡将会引起1、2轨道的积沙均向路堤右侧流动；b-b、c-b路堤的右边坡均设置为下凹坡，左边坡为上凸坡会导致沙粒在1、2轨道沉积，且迎风坡坡脚积沙急剧增加；a-b、c-b路堤的右边坡均设置为下凹坡，左边坡为斜线破或上凸坡均会在1、2轨道出现积沙.由此可知，两侧边坡为下凹坡可减弱风沙运动阻力，有利于减少轨道积沙，即b-b路堤在输沙方面有一定的优势.

为了准确对比各路堤的积沙量，本文在FLUENT中设定路堤顶部1.0 m范围为沙粒体积分数的统计区域，提取、整理数据后，利用Origin软件绘制出四种路堤周围的沙粒体积分数的柱状图(图8).由图8可知，a-b、c-b路堤的沙粒体积分数为0.134、0.147，比a-a路堤的0.120分别高出11.7%和22.5%，表明a-b、c-b路堤的沙粒沉积较多，风沙危害较严重；而b-b路堤的沙粒体积分数最低，比a-a路堤低17.5%，值为0.099，验证了图6(a)轨道积沙较少的结论.综上所述，b-b路堤能够减弱气流运动阻力，减少路堤积沙，降低轨道沙埋的风险，有利于风沙灾害的防护.

图9 四种路堤顶部的沙粒体积分数Fig.9 Sand volume fraction of sand on top of the four embankments

3.3 摩阻风速对路堤积沙的影响

由于风场与运动沙粒的相互耦合作用，风力一定时输沙率也将是特定值，通常用摩阻风速的幂函数来表征输沙率[29]，图10给出了摩阻风速对b-b路堤顶部沙粒体积分数的影响.

由图10可知，当摩阻风速为0.448 m/s时，风的携沙能力较弱，但路堤上回流风速较低，无法再次起动沙粒，致使沙粒沉积较多；当摩阻风速为0.668 m/s、0.892 m/s时，风的携沙能力中等，然而路堤上出现较强回流，风速较大，带动沙粒再次起动远离轨道；当摩阻风速为1.116 m/s时，风的携沙能力较强，同时路堤出现多个较弱回流，轨道附近风速较低，已沉积沙粒无法再次起动，导致积沙较多.综上，随着摩阻风速的增加，b-b路堤顶部的沙粒体积分数先减小后增大，关键原因在于路堤顶部轨道处是否形成了能够带动沙粒再次起动的较强回流.因此，b-b路堤在中等风速范围有明显的输沙优势.

4 结论

沙漠铁路必然会遭受风沙危害，为了保障列车安全运行，必须在路堤两侧增设防风沙措施.目前研究主要集中在挡风墙、挡沙墙等风沙防护措施的防护效益方面，针对路堤自身结构风沙防护效果的研究较少，因此缓和式路堤的风沙流特性研究尤为重要.本文基于气固两相流理论，对比分析了四种路堤周围的流场分布和沙粒沉积状态.

1) 四种路堤周围均形成减速区、加速区及回流区；距路堤0.1 m高度处b-b路堤的最小风速是a-a路堤的2.3倍，值为-9.1 m/s，气流在此处形成较强涡流，风速已达到沙粒起动风速；b-b路堤在路堤背风侧E、F之间区域形成较大回流区，最小风速值为-11 m/s.

2) 与其他模型相比，b-b路堤顶部的积沙最少，比a-a路堤低17.5%；两侧下凹坡的b-b路堤有利于沙粒向远处输运；当摩阻风速为0.668 m/s、0.892 m/s时，b-b路堤的路堤上形成了较强回流，带动已沉积沙粒再次起动远离轨道，降低了沙粒体积分数.

3) 与普通路堤相比，两侧下凹坡的缓和式路堤能加快风沙流通过路堤，降低轨道积沙，是一种输沙性能良好的路堤结构.此外，下凹坡路堤在一定程度上可以防止牲畜穿越轨道影响列车运行安全.新修路堤在考虑轨道积沙量少和实际风速为中等的情况下，建议选择两侧下凹坡的缓和式路堤.