沙漠化不仅威胁着环境，也阻碍着全球的经济发展^[1]，中国是世界上沙漠化最严重的国家之一^[2]；沙漠面积仅小于澳大利亚和沙特，主要分布在西北、华北和东北的部分地区，包括塔克拉玛干沙漠、腾格里沙漠等12个沙漠群，面积约占国土总面积的8.4%^[3]。中国西北地区资源丰富，为推动区域发展，乌玛(乌海-玛沁)高速公路腾格里沙漠段的修建被提上日程。在沙漠地区线路设计时，经常受到地形、不良地质、重要建筑物的影响，须设置平曲线绕避，而曲面迎风形式和曲线半径的不同对风沙活动具有重要的影响，为此诸多学者从不同角度进行了大量研究工作。

近年来，众多学者^[4-7]从行车安全角度出发，利用试验、数据采集和建模等手段系统研究了曲线半径与行车速度、驾驶员心率变化之间的关系，得出了戈壁荒漠地区公路曲线半径的最小取值及推荐取值。不同时期的设计工作者也通过资料总结得出沙漠地区铁路选线的基本原则，净文常^[8]以神榆铁路为研究对象，详细说明了沙漠地区铁路选线工作中应充分考虑风沙活动的影响；李响^[9]总结分析包兰铁路选线资料，发现在沙漠地区选线应尽可能绕避流动沙丘，且平曲线易采用大半径曲线；冯德泉^[10]总结库尔勒-格尔木铁路的选线经验，提出了风沙地区选线应首选曲线外侧朝主风向的建议；马世杰^[11]在和田-若羌铁路选线中提出以桥带路方案通过流动沙丘地区。此外，姚建斐^[12]通过调查分析影响公路沙埋的主要因素，认为大半径曲线由于具有更为舒缓的线型，因此大半径凸面迎风平曲线路堤对风沙流具有较强的疏导作用。

以上学者虽然详细阐述了沙漠戈壁地区平曲线半径的取值范围和铁路选线的基本原则。但平曲线路堤与风沙流之间的作用机理需进一步研究。因此，国内外学者^[13-14]采用数值模拟的手段来研究风沙流与路堤之间的作用机理。2008年，武生智等^[15]首次采用Fluent欧拉双流体模型阐述风沙流与路堤之间的作用机理，并对风沙流作用下的路堤参数做出优化；李驰等^[16]利用风洞试验和数值模拟研究了风沙流作用下路堤、路堑的流场分布特征；石龙等^[17]在此基础上，采用欧拉双流体模型分析了风沙流经路堤时的流场变化及积沙特征。目前，风沙流与路堤的响应规律研究主要以风洞试验和数值模拟为主，且主要从二维角度进行分析，但风沙流结合平曲线路堤的三维数值模拟还少有涉及。

鉴于此，为更好地阐述平曲线路堤与风沙流之间的响应规律，本研究以腾格里沙漠腹地为研究背景，乌玛高速试验段公路路堤为研究对象，基于Fluent欧拉双流体模型，对不同曲面迎风形式、不同曲线半径路堤进行数值模拟，系统地研究了曲面迎风形式和曲线半径对风沙流场结构的影响，结合试验段公路路堤两侧积沙进行验证，进而得出沙漠地区曲面迎风形式和曲线半径的推荐范围，为沙漠地区公路选线提供参考与借鉴。

根据乌玛高速试验段公路路堤横断面参数建立三维简化模型，路堤顶面宽度13 m，路堤高度H=4 m，边坡坡率为1∶2，计算曲线半径分别采用400，1 000，3 000，5 000 m和无穷大(即直线)，曲面迎风形式为凸面迎风和凹面迎风平曲线，线路中心点路堤坡脚记为模型原点，模型如图 1所示。

图 1 模型示意图(单位：m) Fig. 1 Schematic diagram of model (unit: m)

图选项

计算域为x×y×z为300 m×280 m×30 m，路堤模型距离入口60 m，路堤模型距出口最小长度为160 m，路基高度为30 m；参照已有文献资料[18-20]，为保证流场充分发展，模型距计算域上边界高度最小为5H(H为路堤高度，下同)，距入口、出口分别大于10H，20H。因此计算域尺寸合理，对仿真结果无影响。

采用网格划分软件mesh对计算域进行网格划分，网格形式采用四面体与六面体混合形式(如图 2所示)，网格数量大于2 000万，网格质量良好，满足计算要求。

图 2 网格划分示意图 Fig. 2 Schematic diagram of mesh division

图选项

介质类型为流体，截面x=-60 m为速度入口(VELOCITY-INLET)，截面x=220 m为自由出口(OUT-FLOW)，截面z=30 m，y=150 m与y=-150 m采用对称边界条件(SYMMETRY)，其他面采用壁面条件(WALL)，壁面粗糙度默认为0。

研究区沙粒粒径分布区间为0.075~0.25 mm，设定风沙流中的沙粒粒径为0.1 mm，沙粒密度为2 650 kg/m³，沙粒运动黏度系数^[21]为0.047 kg/ (m·s)，入口风速采用均匀风速，模型入口处风速为15 m/s，沙粒相体积分数为0.02。

目前，风沙流的模拟方法有欧拉-欧拉法(双流体模型)和欧拉-拉格朗日法(离散相)。欧拉双流体模型是将风相与沙相分别考虑为连续介质，充满整个流场，用体积分数描述各相在流场中的分布，各项体积分数之和为1，其优点：(1)计算量相对较低；(2)可以充分考虑沙相自身的湍流扩散和风相与沙相之间的速度差与相互作用，使得模拟的结果更贴合实际的沙漠风沙流流动。欧拉-拉格朗日法是将风相考虑为连续介质，采用欧拉方法进行求解，将沙相考虑为离散介质，采用拉格朗日方法进行求解，得到大量颗粒在流场中的运动轨迹及运动参数，但对计算机性能要求高，计算量较大。

沙漠风沙流一般处于饱和状态，沙粒体积分数较大，宜用双流体模型模拟，而对于戈壁风沙流一般都处于非饱和状态，沙粒体积分数较小，宜用欧拉-拉格朗日法。目前大多数研究者就沙漠风沙流普遍采用欧拉双流体模型^[22]。

因此计算模型采用欧拉双流体模型，求解算法采用simplec算法。研究区风速不大于50 m/s，且沙粒主要受重力影响，风沙流可视为不可压缩流体，风沙流之间热量交换可忽略不计，故不考虑能量方程。

欧拉双流体模型模型视计算域内的流固两相物质为互相贯穿的连续介质，两相体积分数之和为1，气相和沙相分别由各自的质量和动量守恒方程控制。流体在近地面或者地形不规则处会出现运动复杂的涡流，多个涡流叠加起来就会形成湍流，气固两相流满足标准k-ε模型的湍动输运方程。

选取线路中间点y=0 m处截面，研究不同曲面迎风形式、不同曲线半径路堤周围风沙流场活动规律，结果如图 3所示。从图(a)~ (e)中均可看出，风沙流途经路堤时，速度发生明显的分区，风沙流在迎风侧坡脚、坡面、路肩分别形成减速区、加速区、高速区，在背风侧坡脚由两个不同涡流形成紊流区，且存在速度突变；在不同曲线半径的路堤形式下，均存在类似情况，可见曲面迎风形式及半径取值大小不影响路堤横截面速度分区为厘清沿线路不同位置处的平曲线路堤风沙流场活动规律，取凸面迎风平曲线、曲线半径为1 000 m时，路堤y=-140 m，y=0 m，y=140 m的3个截面进行流场分析，结果如图 4所示，发现3处路堤横截面流场分布相近，并没有太大改变，说明半径相同时，路堤各个截面位置处的流场基本相同，路堤横截面流场变化与曲线半径大小及所处不同截面位置无关，仅与路堤自身特点^[15-16]有关。

图 3 不同曲线半径路堤横截面流场 Fig. 3 Cross section flow field of embankment with different curve radii

图选项

图 4 R=1 000 m不同截面处路堤流场 Fig. 4 Flow field of embankment at different sections when R=1 000 m

图选项

在背风侧，对不同半径曲线路堤距地面高度1 m(z=1 m)截面进行流场分析，结果如图 5所示，背风侧形成涡流，风速由0恢复至起沙风速并逐渐增大；背风侧涡流区随着半径的增加而减小。

图 5 不同半径下z =1 m截面处路堤流场分布 Fig. 5 Flow field distribution at embankment section with different radii when z =1 m

图选项

当风速恢复至初始风速15 m/s的40%和60%时，该点距路堤背风侧坡脚的长度记为风速恢复区长度，如表 1所列，风速恢复区长度随半径增大而减小。不同曲面迎风形式下，当曲线半径不小于3 000 m时，风速恢复区长度基本保持不变，分别为46，90 m和44，90 m。

对迎风侧z=1 m处截面，提取x=-3，-5，-10，-15，-20，-30 m的风沙速度，沿线路方向进行分析，结果如图 6所示，对凸面迎风平曲线路堤而言，风速沿线路方向呈现“V”形，线路中间点处风速较小，两侧风速较大，说明凸面迎风平曲线路堤对风沙流有侧向输导作用；对凹面迎风平曲线路堤而言，风速沿线路方向呈现“∩”形，线路中间点风速较大，两侧风速较小，凹面迎风平曲线路堤对风沙流具有汇集作用。

图 6 高度1 m处平曲线水平速度变化 Fig. 6 Horizontal velocity variation on flat curve at height of 1 m

图选项

风速变化率Δv₁反映凸面迎风平曲线路堤对风沙流的侧向疏导作用，风速变化率Δv₂反映凹面迎风平曲线路堤对风沙流的汇集作用。Δv₁越小，输导作用越弱，Δv₂越小，汇集作用越弱。当Δv＜1.15(风速变化率小于直线地段最大风速变化率)时，输导作用和汇集作用均可忽略不计，可认为该长度为路堤对风沙流影响的最远长度。对图 6数据进行处理可得z=1 m截面距线路不同位置处的风速变化率，如表 2所列。

对于凸面迎风平曲线，半径R=400 m时，距离线路3 m(x=-3 m)处速度变化最显著，侧向输导作用最强，速度变化率Δv₁为8.99%；对于凹面迎风平曲线，半径R=1 000 m时，距离线路3 m(x=-3 m)处速度变化最显著，汇集作用最强，速度变化率Δv₂为6.66%。

当半径相同时，距离线路越近，侧向输导作用及汇集作用越强，距离线路越远，路堤对风沙流的影响越弱，侧向输导作用及汇集作用越弱；对于凸面迎风平曲线路堤，风沙流输导作用的影响长度约为15 m；对于凹面迎风平曲线路堤，当半径R≤ 1 000 m时，风沙流汇集作用影响长度约为15 m，当半径R＞1 000 m时，路堤对风沙流影响长度仅在10 m左右。

当观测点与线路距离相同时，随着半径增大，风速变化率基本呈下降趋势，路堤对风沙流输导及汇集作用逐渐减弱。

在迎风侧，凸面迎风平曲线路堤对风沙流有侧向输导作用，凹面迎风平曲线路堤对风沙流有汇集作用；当R≥3 000 m时，凸面迎风平曲线输导作用不明显，当R≥5 000 m时，凹面迎风平曲线汇集作用不明显；而在路堤背风侧，当R≥3 000 m时，风速恢复区长度基本保持不变。

因此，在选线过程中，应首选直线，其次为半径不小于3 000 m的凸面迎风平曲线，最后选择半径不小于5 000 m的凹面迎风平曲线。

在乌玛高速K171+800右侧500 m处设置试验段，利用无人机三维航测，观察不同曲线半径路堤的风沙活动特征，采集该区域气象观测站数据，进行风向、风速分析，结果如图 7所示，冬季起风沙风向以西北风、西西北及北西北为主，合成输沙方向124.99°。凸面迎风平曲线路堤迎风侧按照积沙量不同划分为3个区段，分别为区域1，2，3，可以明显看出，在曲线外侧，两侧积沙量较中间多，说明凸面迎风平曲线路堤对于风沙流具有输导作用；区域4位于路堤背风侧，积沙量较迎风侧多，背风侧形成涡流，风速较低，风沙易堆积。曲线路段现场积沙情况从侧面验证了数值模拟的准确性。

图 7 冬季乌玛高速试验段输沙势 Fig. 7 Drift potential of test section of Wuhai-Maqin expressway in winter

图选项

乌玛高速公路全长122.9 km，在K164+870~K183+000约18 km的线路依次通过了固定沙地、平沙地、流动沙丘、平沙地，极易遭受风沙的侵害。对沿线平曲线半径取值进行调查统计，结果如图 8所示，在乌玛高速沙漠段，最小曲线半径采用2 500 m，最大曲线半径为6 500 m，小半径曲线考虑较少，大多采取较大半径曲线，与本研究计算结果(凸面迎风平曲线半径应不小于3 000 m，凹面迎风平曲线半径应不小于5 000 m)基本一致。

图 8 乌玛高速腾格里沙漠段曲线设计值 Fig. 8 Curve design values of Tengger Desert section of Wuhai-Maqin expressway

图选项

(1) 通过数值模拟可发现，路堤周围风速存在分区现象，路堤横截面流场变化与曲面迎风形式及曲线半径无关。

(2) 凸面、凹面迎风平曲线路堤对风沙流分别具有侧向输导、汇集作用，且随着平曲线半径的增大而减小；现场积沙可侧面验证数值模拟结果的准确性。

(3) 当平曲线半径3 000 m及以上时，路堤背风侧风速恢复区长度基本保持不变。

(4) 沙漠地区选线应充分考虑风沙活动特征影响，平面线形首选直线。其次为凸面迎风平曲线，半径建议不小于3 000 m，最后为凹面迎风平曲线，半径建议不小于5 000 m。