高速公路风吹雪灾害分析和挡雪板参数研究

2023-01-18戴明

北京交通大学学报 2022年6期

戴明

（中国铁建投资集团有限公司，广东珠海 519000）

风吹雪是指风速超过临界起动风速后，由气流挟带起分散的雪粒在近地面向前运动的自然现象，这种典型的两相流运动导致积雪重新分布.风吹雪在全球分布广泛，频繁发生，对自然环境和社会经济的影响很大［1］.道路交通中的风吹雪灾害主要体现在其降低路面能见度和阻碍车辆通行，我国西北和东北等高纬度和高海拔地区为风吹雪现象高发地区，道路交通沿线风吹雪灾害呈现出“点多面广线长”且时空分布不均匀的特点.

国内外关于风吹雪灾害及防治的研究成果主要集中在风吹雪机理等基础理论和交通风雪灾害防治技术应用方面，采用的研究方法有野外观测、风洞试验和数值模拟等.通过野外观测获得了风吹雪发生时的风速风向［2］、温湿度［2］、雪粒浓度［3］、风吹雪频率［4］、吹雪层厚度［5］、雪粒跃移通量［6］、临界起动摩阻风速［6］、风吹雪输运总量和悬移通量与摩阻风速的关系［7］等参数.通过室内风洞试验得出了风吹雪过程中雪粒粒径［8］、风雪跃移运动的临界风速［9］、跃移长度［10］、雪粒与床面的相互影响规律［11］，提炼出风雪跃移运动的函数［12-14］，为风雪跃移运动数值模拟提供了基础.数值模型从经验模型到双流体（E-E）模型［15］，再到风雪跃移运动的颗粒轨迹（E-L）模型［16］.目前基于颗粒追踪的三维风吹雪模型逐步建立，但由于湍流的复杂性，很多工作还需进一步探索.上述研究成果部分应用到公路铁路防雪工程设计中，例如美国怀俄明州高速公路自1971年设置防雪栅后，道路区域积雪减少近75%，每年减少由于积雪导致的道路关闭时间9 d，且显著减轻道路灾害事故量［17］.我国的风吹雪研究在风雪流形成机理、运动规律、时空分布及其灾害防治等方面取得较多进展［18］.进入21世纪以来，高寒地区大量公路铁路等交通基础设施建设导致面临的风吹雪问题逐渐增多，在风吹雪基础理论研究及工程应用研究上发展较快［19-21］.其中，防雪栅栏和挡雪板等措施在新疆艾肯达坂公路、精伊霍铁路、阿富准铁路、克塔公路等交通设施沿线发挥作用.但由于防雪设施的挡雪效果与其高度、底部离地间隙、孔隙率等设计参数有关，不同参数的组合影响会产生不同的积雪沉积效果.

风吹雪的发生会受到当地地形地貌、风速风向、降雪程度等局部气候环境以及路基结构形式等多种因素的共同影响，具有很强的地域性，因此对于已有的相关理论、技术和措施不能照搬或直接使用，需要根据实际情况开展分析研究.

新建京新高速公路（G7）新疆巴里坤至木垒段（以下简称巴木段）位于东天山北麓，长约260 km，按双向四车道高速公路标准建设，设计速度为120 km/h.项目所经部分地段降雪期无人居住和行走，缺乏气象观测资料，基于巴木段所经区域风雪时空分布复杂的现状，本文采用现场监测及数值模拟分析，获取了该路线易发生风吹雪灾害区域内的风雪场和环境条件，总结了不同路基结构形式下的风吹雪规律，进而依据实测雪粒参数建立多相流数值模型，分析了挡雪板结构参数对路基内外风雪场的影响，并给出挡雪板设置参数，为工程建设和运营中风吹雪防护提供技术依据.

1 风吹雪特征

1.1 研究背景与监测条件

京新高速新疆巴木段属于大陆性寒冷干旱气候区，冬季严寒，降雪量大，乱风现象明显，气候条件复杂.根据历史气象和雪害调查资料，巴里坤县降雪最大积雪深度达38 cm，木垒县最大积雪深度为44 cm，整个冬季降雪日数平均在40 d，降雪量94.4 mm.巴里坤和木垒地区风向多为南风、西南风和西风，巴木段路线大致为东西走向，风向与路线交角较大，极易发生风吹雪害.此外，高速公路沿线地形地貌的较大差异也会导致风吹雪害严重程度不一.

综合考虑高速公路沿线的地形、地貌、主导风向、风吹雪危害程度、填挖高度及桥涵工程等因素，重点选择浅路堑、深路堑、互通、路堤路堑过渡段、下穿、隧道等风吹雪较严重路段具代表性的工点，设置22处监测点.监测点位如图1所示.

图1 巴木段监测点位示意图Fig.1 Schematic diagram of monitoring points in Bamur Railway Station

采用自主设计的自动气象站获取监测点位的温湿度、1.5 m和3 m梯度风速以及风向数据，每套监测设备可分为供电保障系统、数据采集系统、数据传输系统和物联网云平台系统等部分.监测时间从2020年11月11日至2021年3月20日，总计130 d.采集数据实时上传至新疆京新高速风吹雪监测物联网平台.监测设备如图2所示.

图2 现场监测设备Fig.2 Field monitoring equipment

1.2 现场风吹雪特征分析

1.2.1 环境条件

环境条件主要包括温度、湿度、光照、地形地貌等条件，温湿度与光照条件会明显改变积雪密度和粒径，粒径直接影响雪粒起动风速，而雪粒会在蠕移、跃移和悬移运动中相互碰撞，磨去其晶体枝桠结构，在数值研究中通常简化为粒径符合正态分布的球体结构［22］.

图3巴里坤县2020年11月至2021年2月气温曲线图Fig.3 Temperature graph of Barkol from November 2020 to February 2021

图3 为巴里坤县气象站2020年11月至2021年2月气温变化数据，期间平均气温在-15℃～-10℃之间，木垒县气象站与巴里坤气象站温度变化趋势一致.文献［23］指出雪粒起动速度与环境温度呈指数关系，在温度低于-6℃时雪粒起动风速约为3.4 m/s.沿线区域气候干燥，导致低温条件下雪粒表面难以冻结形成冰壳，在风驱使下易处于运动状态.在巴木段沿线取4个测点测量雪后积雪密度、雪粒起动风速、雪粒粒径等后续试验必要参数.其中，积雪密度采用称重法，所得结果如表1所示.可见高速公路所在区域沉积雪与我国其他地区冬季实测积雪密度接近，新雪密度则明显小于其他地区，因此后续数值计算雪粒密度选定为4个测点的平均值0.127 g/cm3.采用可移动风洞对现场风吹雪的起动风速进行测定，当高度为5 cm测点（对应实际高度3 m）的风速为3.37 m/s时，雪粒产生明显运动，考虑到目视误差，后续数值计算雪粒起动阈值风速设定为3 m/s，如图4所示.

图4 风速测量设备和风速测量剖面Fig.4 Wind speed measuring equipment and wind speed measuring profile

表1 积雪密度测量结果Tab.1 Measuring results of snow density

1.2.2 风场条件

风是风吹雪现象的动力来源，以3 m/s为雪粒起动风速，K219+540监测点降雪期间的30 min平均风速大于起动风速的占比46.9%，图5为该监测点起动风速玫瑰风图.由图5可知，该区域盛行西北风和北风，分别占起动风速的52.3%和36.2%，路线呈东-西走向，即区域内主要起动风向与路线呈垂直或斜交；监测点最大瞬时风速在8.3 m/s，平均风速2.55 m/s.从风速和风向角度均可以看出该监测点位区域易出现风吹雪现象.

图5 K219+540监测点大于起雪风速玫瑰风图Fig.5 Rose wind map of wind greater than starting snow speed during snowfall at K219+540

1.2.3 雪场条件

雪是风吹雪的物质来源，2021年1月19日，巴里坤县和木垒县降水量分别为10.5 mm和14.5 mm，雪后对巴木段里程K198+700～K259+806进行了雪情调查.

K216+000附近路面积雪厚度在50 cm以上，路堤左侧为来流方向，如图6所示.高速公路波形梁护栏和中央隔离带对风雪流运动产生较大影响，无路侧护栏或缆索护栏路段的风雪流通过阻碍小，路面近乎无积雪.

图6 K216+000附近路段积雪情况Fig.6 Snow accumulation along K216+000 section

2 风吹雪与路基形式关系

2.1 路堤试验段

连续填方路段K206+100断面的高度为5 m，路堤-路堑过渡段K219+540断面的路堤高度为1 m，分别测试其风速流场和积雪分布，结果如图7和图8所示.由图7可以看出K206+100断面上风侧边坡坡脚处风速最小，沿路堤边坡向上风速逐渐增高，至路肩附近达到最大.路面风速大于来流风速.下风侧边坡风速顺坡向大幅度下降，至坡脚处风速降低到最小值.路基边坡积雪面倾角为56°，符合雪粒子休止角的范围.路堤断面各处风速分布的不均匀性直接影响到积雪分布特征，雪粒主要堆积于路堤两侧坡脚位置，如图8所示，积雪深度最大达到71 cm，路面雪深仅为2～3 cm.

图7 路堤断面风速分布Fig.7 Distribution of wind speed at the cross section of embankment

图8 路堤积雪雪深分布Fig.8 Distribution of snow depth at the cross section of embankment

K219+540断面所在区域被积雪完全覆盖，道路轮廓不清楚，受风吹雪影响很大.路堤各部位风速变化情况与K206+100断面一致，但路面风速较小且波动幅度不明显.下风侧坡脚积雪厚度达到101 cm，路面积雪平均厚度为12 cm，分别是K206+100路段的1.4倍和4.8倍.由此可见，风速对积雪分布的影响显著.风速越低，产生的驱动力不足以平衡雪粒之间和雪粒与地面之间的摩阻力，往往产生较为明显的积雪分布，因此工程中应注意对低风速地区的防雪除雪.此外，低路堤易受到路侧障碍物的影响，应清除路基两侧凸出的乱石堆、小土丘、草墩、灌木丛等障碍物，避免路基或路面范围出现低风区.

2.2 路堑试验段

分别测试填挖过渡路段K164+000断面和连续路堑路段K219+940断面的风速流场和积雪分布，结果如图9和图10所示.K164+000断面路堑深度约5 m，从上风侧边坡开始风速下降，在坡顶附面层分离，风速急剧减弱.上风侧路肩处风速分别为1.4 m/s和1.7 m/s，仅为来流风速的46%和32%，雪粒子堆积在上风侧边坡，积雪形成沉积锋面，在锋面前缘形成新的减速区，并逐渐向路面推进.位于上风侧边坡坡脚和路基坡脚之间的左侧积雪平台处风速变化明显，积雪较厚，最大积雪深度为35.5 cm.路面平均风速分别为1.5 m/s和2 m/s，明显低于来流风速，风雪流经过路面时有少量雪粒子沉积.下风侧边坡风速再次减小，在右侧积雪平台处形成减速区，该处最大积雪深度为24 cm.风雪流到达坡顶后，风速逐渐恢复.

图9 路堑断面风速分布Fig.9 Comparison of wind speed at the cross section of cutting

图10 路堑积雪雪深分布Fig.10 Distribution of snow depth at the cross section of cutting

K219+940与K164+000断面路堑深度接近，但是前者路堑内积雪更为严重，分析原因是其位于长度300 m的连续路堑，地形对风雪流的扰动范围远大于连续填挖过渡路段路堑.

3 数值模型的建立与验证

3.1 运动模型与方程

本文采用不可压缩流雷诺时均方程（Reynolds-Average Navier-Stokes，RANS）表示流场特征，其动量方程为

在其中加入重整化群的k-ε湍流模型，以封闭RANS方程中的时均项进行求解，其中湍流模型的湍动能k为

式中：t为计算时间；ρ为空气密度，取值为1.225 g/cm3；p为压强；xi和xj表示流场和重力方向上的长度分量；ui和uj为平均速度分量；u′i和u′j为脉动速度分量；σij为应力张量分量；μ表示空气动力黏度，取值为17.9×10-6Pa/s，μt由Cμ计算得出，Cμ取值为0.084 5；Gk表示层流速度梯度产生的湍流动能；αk和αε分别表示ε方程的湍流普朗特常数和湍流额外波动，取值均为1.39；C1ε和C2ε分别为经验系数，取值为1.42和1.68.

上述方法对小范围内流场剧烈变化的情况有较好的模拟精确度［24-25］.

将雪粒视为连续体，采用欧拉多相流模型，引入雪粒体积分数的概念，添加雪相运动控制方程如下，

式中：ρs为雪粒密度，设为现场实测值0.127 g/cm3；fs为雪粒体积分数；Sc为施密特数，取值为1.

公式中的每一项表示雪粒的不同运动状态，其中等式左边第一项为非定常项，表示与时间相关的物理变化量，在定常计算中可忽略该项，等式左边第二项和右边第一项共同表示雪粒运动的对流-扩散特征.计算中各相体积分数之和恒为1，将风吹雪视为不可压缩流，采用Syamlal-Obrien模型考虑气体对雪粒的曳力.

3.2 计算域模型建立

建立高速公路路基与挡雪板模型如图11所示，路基面宽度27 m，边坡坡度1∶4，挡雪板具体尺寸随工况变化，不同工况挡雪板尺寸如表2所示.为保证进口流场的充分发展和减少计算域边界壁面对主要模型域的影响，计算域整体尺寸较主要模型域至少扩大10倍，设定流场方向1 000 m，线路方向500 m，高度50 m.

表2 不同形式挡雪板计算工况Tab.2 Calculation conditions for different forms of snow guards

图11 路基与挡雪板区域模型Fig.11 Mixed mesh of snow snow guard

挡雪板处采用非结构化网格，其他区域采用结构化网格，在保证网格质量的同时减小网格数量，提高计算效率.为考虑极端风雪灾害，假定来流风向均垂直于挡雪板与线路方向.计算模型的边界条件如表3所示.

表3 边界条件Tab.3 Boundary conditions

式中：v(z)为高度z处的风速；u*t为阈值摩擦速度，取值为0.26 m/s，对应流场入口高度3 m处风速为6 m/s；κ为卡曼常数，取值为0.4；z0为壁面粗糙度高度，取值为3.0×10-4m.

计算中采用SIMPLEC算法求解不可压缩流，实现耦合代数方程组的解耦及速度场和压力场的修正，其中压力差值采用Standard格式，其他选用精度较高的二阶迎风格式.

3.3 结果可靠性检验

3.3.1 网格无关性分析

本文采用的Fluent 18.0软件中求解质量方程和动量方程均是以有限体积法为基础，其单元节点位于网格边的中点及网格体的中心，因此模型网格单元尺寸与数量会直接影响结果.由于风吹雪现象涉及明显的湍流作用，因此近地面和主要模型近壁面处采用边界层网格，壁面函数计算要求第一层网格尺寸y+在30～300之间，本文模型y+为34.05，符合计算要求.

表4为不同网格尺寸下的计算条件和计算结果，选取4组不同的网格尺寸，每组中最大结构型网格尺寸均设置为最大非结构型网格的10倍，各组模拟时长均为300 s以保证流场的充分发展和积雪的稳定分布.可以看出当稀疏处结构型网格尺寸为0.5和0.2时，相邻网格间误差值最大达到18.0%和9.6%，当稀疏处结构型网格尺寸为0.1和0.05时，误差值低于5%，4组网格尺寸中加密处非结构型相邻网格间的误差值均小于3%.综合结构型和非结构型相邻网格的计算误差，可以认为第3组和第4组网格尺寸的计算结果精度已满足要求，且随着网格数量的增加，计算

表4 不同网格尺寸下计算条件与结果Tab.4 Calculation conditions and results under different mesh sizes

结果与网格不再保持显著的敏感性.

3.3.2 计算时间相关性分析

在多相流瞬态计算中积雪分布与计算时间呈现出相关性，图12为不同计算时间下无孔隙3m高挡雪墙两侧的积雪分布.对应图12中的雪粒体积分数数据见表5.由表5可知，计算初期，挡雪板两侧雪量堆积迅速增加，随着模拟时长的增加挡雪板两侧雪量也会不断增长，但增长率越来越慢，同时出口边界的雪量则迅速增加.综合考虑积雪分布稳定性和计算时长条件，本文选择t=300 s作为后续模拟时长.

表5 各点位雪粒体积分数Tab.5 Volume fraction of snow particle at different points

图12 不同模拟时长下的雪量分布Fig.12 Distribution of snow under different simulation durations

为了提高计算效率，计算中首先开展单向流场模拟，待其达到收敛后，添加雪粒进行多相流瞬态计算.瞬态计算中时间时长为5×10-2s，单时间步长内迭代次数为20，保证每一个时间步长内残差曲线均达到收敛.

根据上述试算与分析，综合考虑计算时长，最终选定计算域稀疏处结构型网格最大值为0.1，加密处非结构网格最大值为0.01，网格间取1.05的递增率进行加密，模拟时长为300 s.

3.3.3 路堤段对比分析

将K206+100路堤断面风速和雪深的数值计算结果与实际值对比，如图13和图14所示.可以看出，二者变化趋势较为一致，路堤路基风场廓线形状近似于“n”形，积雪分布则与风速分布相反呈“u”形分布.由于路堤的阻碍作用，路基上风侧坡脚处风速剧烈衰减，顺着坡面风速持续增加至坡顶，风速到达路面中央略有减小，道路表面积雪少.背风坡坡脚因存在涡旋减速区，风速再一次衰减，路基两侧雪深较大.

图13 K206+100路堤断面实测风速与数值计算结果对比Fig.13 Wind distribution comparison of measured value and numerical results at K206+100 embankment cross section

图14 K206+100路堤断面实测积雪深度分布与数值计算结果对比Fig.14 Snow depth distribution comparison of measured value and numerical results at K206+100 embankment cross section

3.3.4 路堑段对比分析

对比分析K164+000路堑段风速和雪深的数值计算结果与实际测量值，如图15和图16所示.可以看出两者变化趋势较为一致，路堑路基风场廓线形状近似于“w”形，路堑内积雪分布表现为“m”形.顺来流方向，路堑边坡坡顶风速增大，然后在积雪平台内风速衰减，而道路表面的风速相对积雪平台则有所增大，这种风速变化趋势使得路堑外吹向路堑内的雪粒较多地沉积于积雪平台，有效地减少了道路表面的积雪量.

图15 K164+000路堑断面实测风速与数值计算结果对比Fig.15 Wind distribution comparison of measured value and numerical results at K164+000 cutting cross section

图16 K164+000路堑断面实测积雪深度分布与数值计算结果对比Fig.16 Snow depth distribution comparison of measured value and numerical results at K164+000 cutting cross section

4 挡雪板设计参数分析

4.1 挡雪板底部间隙对风雪场影响

为避免挡雪板底部间隙d的变化引起整体孔隙率的变化，模型中挡雪板孔隙率β设为0，挡雪板高度H为4 m.图17为挡雪板不同底部间隙工况下的风速流场分布云图.可以发现，在d=0时，挡雪板上风侧流场受到完全的阻碍作用，底部速度几乎衰减至0；下风侧出现单个尺度相对较大的涡旋，其中心风速减弱最为明显，而在涡旋外侧风速恢复较快.随着d增大，气流的压缩作用使得近地面存在流场加速区，挡雪板下部流速明显增强，上风侧雪粒更多地通过挡雪板进入下风侧，同时形成的涡旋数量逐渐增多，分布范围也越来越广泛.

图17不同挡雪板底部间隙下风速云图Fig.17 Contour of wind speed under different bottom gaps of snow guard

图18 为挡雪板不同底部间隙工况下的积雪分布云图.可以发现，随着底部间隙从无到有并逐步增大，积雪逐渐由挡雪板前向挡雪板后移动，并且积雪堆积趋势越来越平缓.当d=0时，越过挡雪板上方的雪粒较少，积雪主要集中在挡雪板前大约8m左右位置；随着d增大，挡雪板后并紧靠挡雪板的位置积雪相对无底部间隙时减少，原因是涡旋影响积雪的分布.

图18挡雪板不同底部间隙下雪量体积分数云图Fig.18 Contour of snow volume fraction under different bottom gaps of snow guard

4.2 挡雪板高度对风雪场影响

图19为挡雪板不同高度下雪粒体积分数云图，其中β=0.25，d=0.2 m.可以发现，不同高度挡雪板附近雪粒体积分数最大值随挡雪板高度的增加而向板后偏移，后移距离约为2H.

图19 不同挡雪板高度下雪粒体积分数云图Fig.19 Contour of snow volume fraction under different heights of snow guard

图20为不同高度挡雪板工况下的水平风速分布.在低于挡雪板高度（距离地面0.2 m高）的区域，板后减速效果均较为明显，并且挡雪板后的涡旋区范围随挡雪板高度的增加而增大.随着挡雪板高度的增加，涡旋中心逐渐右移，积雪覆盖范围逐渐增大.H=2 m时，挡雪板的下风侧存在长约25 m，高1 m的涡旋区范围；H=3 m时，挡雪板的下风侧存在长约30 m，高约2.5 m的涡旋区；H=4 m时，挡雪板下风侧存在长约40 m，高约3 m的涡旋区.当监测位置高于挡雪板高度时（距离地面4 m高），板后减速效果随挡雪板高度的增加而愈加明显，在高于挡雪板H时，板后风速变化较小，即挡雪板阻风的影响高度约为H.

图20 不同高度挡雪板两侧流速分布Fig.20 Velocity distribution on both sides of snow guard with different heights

不同高度挡雪板工况下雪量分布如图21所示.可以发现，挡雪板越高，有效阻雪距离及阻雪量越大.β=0.25时，不同高度挡雪板下风侧积雪沉积范围在10H～14H之间.H=2 m时，雪粒体积分数在20～30 m位置逐渐稳定，H=3 m和4 m时，雪颗粒体积分数分别在30～40 m、40～50 m位置减缓.

图21 挡雪板高度对雪量的影响Fig.21 Effect of snow guard height on snow volume

4.3 挡雪板孔隙率对风雪场影响

孔隙率为挡雪板区别于实心挡雪板（β=0）最主要的结构特征.图22为不同孔隙率挡雪板的雪粒体积分数云图，其中H=4 m，d=0.2 m.可以发现由于低孔隙率挡雪板附近风场减速效果明显，雪粒沉积往往较多.β=0的挡雪板的积雪主要分布在上风侧以及下风侧距离挡雪板较近的位置，积雪高度存在一个明显的峰值，但积雪距离较短.β=0.25的挡雪板上风侧的雪量有所降低，下风侧积雪长度有所增加.随着孔隙率的增加，积雪更多的向下风侧堆积，积雪分布范围更为广泛，积雪形态也更加平缓.

图22 不同挡雪板孔隙率下雪粒体积分数云图Fig.22 Contour of snow volume fraction under different porosities of snow guard

不同孔隙率挡雪板风速分布如图23所示.可以发现，实心挡雪板下风侧对风场的影响最为显著.在挡雪板下风侧存在一个长约35 m，高约3 m的风速衰减区，其中最大回流速度为2.6 m/s，在40 m之后风速快速回升；当β=0.25、0.4和0.6时，挡雪板下风侧风速衰减相对减缓，但风速回升点位达到了挡雪板下风侧50、60和80 m左右，即随着β的增大，涡流由小范围的剧烈作用转变为大范围内的平缓影响，挡雪板下风侧的涡旋区中心向更远处偏移，减速范围增大，最大回流速度减小.

图23 挡雪板孔隙率对风速的影响Fig.23 Effect of snow guard porosity on wind velocity

不同孔隙率挡雪板工况下雪量分布如图24所示.可以发现，实心挡雪板的下风侧积雪沉积范围在35 m左右，随后雪粒体积分数趋于平缓；β=0.25和0.4时板后积雪沉积范围约为45 m（11H）、60 m（15H），这与β=0.25时挡雪板高度H和风雪分布特征研究的结论一致.即随着β增加，板后积雪分布曲线变化波动减缓，雪粒沉积范围增大；挡雪板下风侧需要增加储雪距离，避免路基面与挡雪板过近使得路基面处于风速减弱区.