青藏铁路东起西宁, 西至拉萨, 全长1 956 km, 是世界海拔最高的高原铁路, 途经河谷、退化草原、戈壁等多种地貌类型(谢胜波等, 2014), 穿越海拔4 000 m以上地段达960 km, 最高点海拔5 072 m。青藏高原地理位置特殊, 自然灾害频发, 其中风沙灾害是影响铁路安全运行的最重要因素之一(邱青云等, 2004)。青藏铁路多处经过风沙地段, 沙害面积约880 km², 风沙段长度总计442.6 km, 占线路总长的22.5%(Zhang et al., 2010; 张克存等, 2011; 薛智德等, 2010), 其中轻度沙害路段282 km, 中度沙害路段105 km, 严重沙害路段53 km(于云江等, 2001; 张克存等, 2010)。目前, 在沙害严重路段已设立系统的防沙措施并取得良好成效, 但随着时间的推移, 部分防沙措施已逐渐失去作用, 沙害问题仍十分严峻(牛清河等, 2009)。

关于风沙流的研究手段主要有2种, 一种是传统的野外观测及风洞试验, 另一种是近年来兴起的数值模拟。与传统手段相比, 数值模拟不仅可以减少野外的工作量, 而且节约成本, 又可以快速得到结果, 因此越来越多地应用于风沙运动的研究。李驰等(2010)利用Fluent软件对风沙环境下铁路路基的风蚀破坏规律进行了数值模拟研究, 分析了路基高度和边坡比对流场的影响, 结果显示风速减弱区和恢复区主要受路基高度影响。石龙等(2014)基于Fluent软件, 对兰新铁路路堤周围风沙两相流运动特性进行了模拟研究, 揭示了迎风坡积沙量大于背风坡的形成机制。李晓军等(2016)对兰新铁路既有挡风墙周围风沙两相流运动特性进行了数值模拟, 阐述了挡风墙背风侧风速廓线变化规律, 提出了挡风墙建设高度的合理建议。然而, 对于青藏铁路路基流场及积沙的研究仍然十分薄弱, 路基周围流场及积沙分布规律仍然不很清楚。本文通过Fluent软件进行模拟计算, 研究青藏铁路原有路基及新提出的一种输沙型路基的风场分布情况以及路基积沙特征, 并讨论了轨枕式挡沙墙与PE (Polyethylene)阻沙网对路基风场及积沙特征的影响, 探索路基对风沙运动规律的影响, 为青藏铁路防沙治沙提供理论基础。

1 模拟方法

数值模拟方法包括前处理、计算、后处理3个步骤。前处理使用Gambit 2.4建立几何模型及网格划分, 计算过程采用Fluent 16.0, 后处理使用Tecplot 360。具体步骤见图 1。

1.1 几何模型

风沙流是风成沙在多种作用力下的三维运动现象, 需要考虑多个因素在3个空间方向上的变化, 十分复杂, 由于风沙流在运动过程中主要受到水平和垂直方向上力的影响, 因此采用二维模型进行模拟。

本研究对3种不同的路基进行模拟。第1种是普通路基, 根据《铁路路基设计规范》(TB10001—2005), 按比例对普通长直型铁路路基进行建模, 具体规格见图 2a。第2种是通风路基, 根据刘志强等(2006)的研究, 通风路基规格如图 2b, 其中区域Ⅰ为道碴层, 区域Ⅱ为碎石砂砾。通风管在区域Ⅱ中, 管中心线距地面1 m, 管直径0.4 m。第3种是输沙型路基, 根据现有普通长直路基创新而来, 主要适用于主风向单一并与铁路线夹角较大或垂直的严重风沙路段。具体规格如图 2c所示, 其中：进口段的底面斜坡坡度为1/8, 集风段底面斜坡坡度为1/15, 输沙段底面为水平面, 出口段的底面斜坡坡度为1/10, 集沙池在沿路基延伸方向上的底面斜坡坡度依地势分段设置坡降, 其长度约为输沙洞长度的1/4~2/5。

1.2 边界条件设置

计算域左侧边界是风沙流入口, 定义为速度进口边界条件(VELOCITY_INLET), 右侧是风沙流出口边界, 定义为完全发展出流边界条件(OUT_FLOW), 计算区域顶部采用对称边界条件(SYMMETRY), 路基模型及地面定义为无滑移壁面(WALL), 介质类型为FLUID。

本模拟中入口速度通过udf函数自定义, 垂直风速廓线表示为:

$ V = {V_1}{\left({z/{z_1}} \right)^\alpha } 。$

(1)

式中:z代表任意高度; z₁取值为2, 代表 2 m高度; V₁代表z₁高度处的风速; V代表z高度处的风速; α代表粗糙度指数, 取值0.2。

1.3 参数设置

本文对风沙流进行非稳态模拟, 采用欧拉多相流模型, 由于风沙流经过路基时会产生涡旋流, 假定湍流完全发展, 使用标准k-ε湍流模型, 经计算得湍流强度I=0.02。

空气密度ρ_g=1.225 kg·m^-3, 沙粒密度ρ_s= 2 650 kg·m^-3, 粒径分布主要在0.05~0.25 mm之间, 取0.1 mm, 沙粒相体积分数取0.02。空间离散化都使用一阶迎风格式, 时间步长取0.001 s, 流场求解采用SIMPLE算法。

1.4 基本方程

1) 质量守恒方程  对于不可压缩的流体流动, 密度为常数, 表达式如下:

$ \frac{\partial }{{\partial t}}\left({{\alpha _q}{\rho _q}} \right) + \nabla \cdot \left({{\alpha _q}{\rho _q}{{\vec v}_q}} \right) = 0 。$

(2)

式中:α_q表示第q相体积分数, 风沙流是气固2相流, 因此气体相与沙粒相体积分数之和为1, 即α_g+α_s=1;ρ_q为第q相的密度; ${{{\vec v}_q}}$为第q相的矢量速度。

2) 动量守恒方程  气体相:

$ \begin{array}{l} \frac{\partial }{{\partial t}}\left({{\alpha _{\rm{g}}}{\rho _{\rm{g}}}{{\vec v}_{\rm{g}}}} \right) + \nabla \cdot \left({{\alpha _{\rm{g}}}{\rho _{\rm{g}}}{{\vec v}_{\rm{g}}}{{\vec v}_{\rm{g}}}} \right) = \\ \;\;\; - {\alpha _{\rm{g}}}\nabla p + \nabla \cdot {\overline{\overline \tau } _{\rm{g}}} + {\alpha _{\rm{g}}}{\rho _{\rm{g}}}\vec g + {f_{{\rm{sg}}}} ；\end{array} $

(3)

沙粒相:

$ \begin{array}{l} \frac{\partial }{{\partial t}}\left({{\alpha _{\rm{s}}}{\rho _{\rm{s}}}{{\vec v}_{\rm{s}}}} \right) + \nabla \cdot \left({{\alpha _{\rm{s}}}{\rho _{\rm{s}}}{{\vec v}_{\rm{s}}}{{\vec v}_{\rm{s}}}} \right) = \\ - {\alpha _{\rm{s}}}\nabla p + \nabla {p_{\rm{s}}} + \nabla \cdot {\overline{\overline \tau } _{\rm{s}}} + {\alpha _{\rm{s}}}{\rho _{\rm{s}}}\vec g + {f_{{\rm{sg}}}} 。\end{array} $

(4)

式中: ${\vec g}$为重力加速度; p为2相共享的相同压力; p_s是沙粒相压力; ${\overline{\overline \tau } _{\rm{g}}}、{\overline{\overline \tau } _{\rm{s}}}$分别为气体相和沙粒相的应力-应变张量; f_sg为气体相和固体相间的相互作用力。

本文使用标准k-ε湍流模型, 由于沙粒相体积分数只有0.02, 因此气体相的湍流运动在风沙流中占主导作用, 假设湍流完全发展, 且不考虑源相的影响。

3) 湍流动能k方程:

$ \rho \frac{{{\rm{D}}k}}{{{\rm{D}}t}} = \frac{\partial }{{\partial {x_i}}}\left[ {\left({{\mu _l} + \frac{{{\mu _t}}}{{{\sigma _k}}}} \right)\frac{{{{\rm{δ }}_\mathit{k}}}}{{{\rm{ δ }}{x_i}}}} \right] + {G_k} + {G_b} - \rho \varepsilon ；$

(5)

4) 湍流耗散率ε方程:

$ \begin{array}{l} \rho \frac{{{\rm{D}}\varepsilon }}{{{\rm{D}}\varepsilon }} = \frac{\partial }{{\partial {x_i}}}\left[ {\left({{\mu _l} + \frac{{{\mu _t}}}{{{\sigma _\varepsilon }}}} \right)\frac{{{{\rm{ δ }}_\varepsilon }}}{{{\rm{ δ }}{x_i}}}} \right] + \\ {C_{1\varepsilon }}\frac{\varepsilon }{k}\left({{G_k} + {G_{3\varepsilon }}{G_b}} \right) - {C_{2\varepsilon }}\rho \frac{{{\varepsilon ^2}}}{k} ；\end{array} $

(6)

式中：μ=μ_t+μ_l, μ_t= $\rho {C_u}\frac{{{k^2}}}{\varepsilon }$。其中：G_k是由层流速度梯度而产生的湍流动能; G_b是由浮力产生的湍流动能; μ是有效黏性系数; μ_l是层流黏性系数; μ_t是湍流黏性系数; C_1ε、C_2ε、G_3ε、σ_k, σ_ε为经验常数, C_u为湍流常数。

2 结果与分析 2.1 路基平面流场速度分布特征

根据观测资料, 青藏铁路部分路段持续时间较长的输沙风速为15 m·s^-1左右, 因此在模拟时分别测试了高度2 m处风速为10、15、20 m·s^-1时的流场情况, 沙粒相体积分数取0.02, 3种路基的流场速度分布见图 3。

由图 3可知, 不同风速的流场没有明显差异。3种路基的流场存在相似分区, 风沙流在平坦戈壁上流动, 途经路基阻碍时, 在迎风侧坡脚处速度降低形成减速区(A); 路基迎风坡在垂直方向上对气流有抬升作用, 速度逐渐增大, 形成抬升区(B); 路基上方由于气流聚集, 风速增大, 形成高速区(C); 在路基背风侧有明显的涡旋流, 风速降低, 部分区域风速变为负值, 该区域风向与来流风向相反, 形成回流区(D); 风沙流远离路基之后, 气流开始扩散, 速度缓慢上升, 最终恢复到初始状态, 为恢复区(E); 此外, 针对通风路基和输沙路基, 由于路基内部存在管道, 气流从开阔地带流入狭窄的管道内, 空气不能大量堆积, 因此加速流过管道, 形成加速区(F)。

2.2 路基边坡特征位置流场速度变化情况

从上述分析中可知, 风速对流场变化的影响较小, 本部分着重分析风沙流速度为15 m·s^-1时, 路基边坡特征位置处的速度变化情况。在野外调查中, 笔者发现积沙主要位于路基迎风侧坡脚、铁轨内和背风侧坡脚, 因此分析这3个断面0~10 m的风速廓线(卢玉华等, 2012)。由图 4可见, 在迎风侧坡脚处, 3种路基的风速廓线变化趋势一致, 输沙路基的风速随高度增幅最大, 说明在迎风坡输沙路基对气流阻碍作用最弱, 沙粒最不易遇阻沉积。路基顶部选取铁轨中轴位置, 在0~2 m范围内, 气流穿过通风路基和输沙路基内部的管道, 风速激增, 管道中部风速高于近管壁处风速, 并且输沙路基的管道内部风速更大; 2~5 m处是路基的实体部分, 风速值为0。路基背风侧存在明显的涡旋流, 流场情况较为复杂, 风速都呈现上升、下降、再上升的趋势; 在0~5 m范围内, 普通路基的风速始终低于5 m·s^-1, 因此最可能导致沙粒在坡脚沉积。

边坡侵蚀和铁轨磨蚀主要是由于近地层气流中沙粒对边坡和铁轨的冲击(邹维, 2012), 因此分析路基特征表面30 cm高度处的风速变化情况。从表 1可以看出, 普通路基4个特征位置处的风速, 都比入口速度(15 m·s^-1)低, 说明在路基表层有一层低速区域, 速度低于起沙风速, 所以普通路基边坡积沙较为严重。通风路基迎风侧坡脚风速减小约64%, 路肩风速增加约19%, 气流在迎风坡上升时, 速度显著增大。输沙路基迎风侧坡脚风速减小约36%, 路肩风速增加约25%, 在3种路基中对气流的阻碍作用最小, 因此不易在路基边坡形成积沙。背风侧坡脚风速减小约20%, 由于输沙路基输沙段促使气流加速, 所以在背风侧坡脚风速仍较高, 经模拟证明该路基在单风向风沙较大的区域具有良好的适用性。

2.3 路基边坡积沙分布特征

青藏高原的风向与铁路的走向有一定夹角, 而且很多地段处于多风向地带, 不同的风向会影响路基的积沙情况。本文主要考虑风向与铁路垂直时(风沙流速度为15 m·s^-1)的坡面积沙情况, 这对多风向的实际情况具有一定的借鉴意义, 并能反映出边坡积沙的普遍规律。

图 5是风沙流经过时路基周围积沙分布图(沙粒的体积分数表示该区域沙粒体积占总体积的比值, 其值越接近1表示该区域积沙量越大。不同颜色代表沙粒不同的体积分数, 红色最大, 蓝色最小, 红色表示沙粒已沉积在此处, 其他颜色代表沙粒仍在继续运动)。从图 5可以看出, 在风沙流经过路基时, 3种路基均有不同程度的积沙, 且路基迎风坡积沙量明显大于背风坡。普通路基边坡积沙最为严重, 迎风侧路堤和路肩都存在积沙现象, 路基近地层气流中沙粒含量较大, 导致风沙流持续吹蚀边坡及铁轨。通风路基的迎风坡脚和背风坡脚均有部分积沙, 通风管上部的边坡没有积沙, 但是通风管内积沙严重, 可能会造成通风管堵塞, 长此以往必定会使其保护冻土的作用减弱, 甚至失效。输沙路基边坡没有积沙, 当风沙流经过路基时, 一些沙粒沉积在输沙段内, 表层沙粒继续运动, 最终在积沙池内沉积, 其余沙粒随气流移动, 离开路基。整个过程中路基边坡以及顶部无沙粒沉积, 可见输沙路基在单风向风沙较严重的地区, 可以起到减轻风蚀的作用, 并且人工清理积沙池内的积沙也相对便利。

2.4 沙障对路基流场及积沙情况的影响

根据青藏铁路南山口路段的实际调查情况, 该路段在路基上风向布设了多种沙障组合, 一共5个分区(姜鑫贵等, 2016)。本部分主要分析轨枕式挡沙墙(以下简称栅栏沙障)和PE阻沙网对路基流场及积沙的影响。

本模拟中栅栏沙障高度2.1 m, 厚度0.3 m, 由废旧轨枕并排埋设而成。PE阻沙网高度为1.5 m, 整个沙障均匀透风, 孔隙度约为40%。从图 6可以看出, PE阻沙网和栅栏沙障布设在上风向时, 风沙流遇阻会有明显的减速过程, 2种沙障对迎风侧风场的影响范围均小于背风侧, 并且PE阻沙网对风场的影响范围略小于栅栏沙障, 这是因为PE阻沙网的高度小于栅栏沙障, 而且PE阻沙网具有一定孔隙度, 气流可以穿过阻沙网(唐玉龙, 2013)。风沙流离开沙障影响范围后, 速度开始缓慢增加, 当抵达路基迎风坡时, 又遇阻下降。与无沙障布设时相比, 风沙流在路基迎风坡的速度明显降低, 但由于沙障将大部分沙粒阻挡, 并且障后风速不足以起动沉积在背风侧的沙粒, 因此路基边坡的积沙量大大降低(杨印海等, 2010)。从模拟结果来看, 2种沙障之间存在一个低速区, 由于阻沙网对风沙流的阻碍作用, 阻沙网背风侧风沙流减速沉降, 随后在平坦开阔的戈壁地表流速有所回升, 但经过栅栏沙障时速度再次下降, 所以这个区域的风速呈现两端小、中间大的趋势, 这个规律在栅栏沙障与路基之间也是存在的, 只是由于两者间距较小, 所以并不明显。

图 7是2种沙障的积沙分布情况, 栅栏沙障迎风侧积沙量大于背风侧, PE阻沙网背风侧积沙量大于迎风侧, 辛国伟等(2016)以及郝晓杰(2012)也发现了这一结论。图 9为南山口沙障积沙情况, 可以看出沙障两侧实际积沙形态与本文模拟基本一致。从沙障的积沙情况来看, 栅栏沙障迎风侧积沙量较大, 并且还有部分沙粒继续以跃移和悬移的形式运动, 但由于风速很低, 最后沙粒会在迎风侧堆积。PE阻沙网迎风侧没有明显沙粒堆积, 大部分沙粒会通过孔隙直接穿越到阻沙网背风侧, 只有少部分沙粒会随气流攀升, 在PE阻沙网上部穿过, 最终沉积在背风侧。对比2种沙障的积沙分布特征可以看出, 栅栏沙障的拦沙效果明显比PE阻沙网好, PE阻沙网背风侧积沙量较少, 堆积高度很低, 有部分沙粒会在风力作用下继续运动。

图 8为布设沙障时路基积沙情况, 与图 5相比, 3种路基边坡的积沙量都相对减少。普通路基迎风坡的积沙情况有很大改观, 只有极少的沙粒堆积在坡脚, 风沙流沿着边坡抬升, 经过铁轨后在背风坡脚沉积, 其余沙粒仍处在运动状态, 随气流搬运远离路基。在通风路基上风向布设沙障后, 有部分沙粒在边坡移动, 但并未堆积, 且通风管内的积沙也明显减少, 贴近管道下部堆积的沙粒层变薄, 保温功能不会受到太大影响。输沙路基是这3种路基中积沙最少的, 整个路基结构几乎没有积沙, 积沙池内沙量较少, 积沙池后也有一小部分积沙, 人工清理相对便利。可以看出沙障的存在限制了沙粒的运动, 大量沙粒被拦截, 只有少部分沙粒被搬运至路基附近, 起到了很好的防沙效果。

3 讨论

本研究模拟结果显示, 随着风速的增加, 3种路基周围流场形态差异不大, 都具有明显的分区。张克存等(2010)对沱沱河路段路基流场结构进行了风洞模拟试验研究, 也发现气流在经过路基时形成明显分区。姚立强等(2013)对临策铁路戈壁段4 m高度路基的流场实测结果与本文模拟的流场分布相近。路基边坡的蚀积特征与气流速度有很好的对应关系, 在路基上方的高速区产生强烈风蚀, 沙粒主要堆积在路基两侧坡脚低速区(张军平等, 2011)。笔者通过模拟发现路基迎风侧坡脚的积沙量大于背风侧, 而不同种类的路基边坡积沙情况也有所不同。尤全刚等(2011)也指出路基背风坡坡度越小越易造成沙粒沉积, 并且随着风速增加沙粒沉积位置逐渐从背风坡中部向坡脚移动。风沙路段路基边坡的损害主要是风蚀破坏, 路基高度在一定程度上影响风蚀以及沙粒的堆积位置(苑宝华等, 2012; 刘伟伟等, 2006)。笔者提出的输沙路基主要起输导作用, 气流经过输沙管道时先迅速增大再缓慢减小, 管道的长宽比会影响风速增幅(洪新民等, 2017), 可根据实际情况适当调整以保证风速大于起沙风, 携沙气流始终从输沙管通过, 因此风沙流对路基边坡的风蚀及对铁轨的磨蚀作用有所减轻, 边坡积沙量明显减少。

高立式沙障的防风固沙作用主要体现在降低风速和促使沙粒从空气中分离并沉积(刘世海等, 2010)。阻沙栅栏是青藏铁路沿线使用最多的防护措施之一, 主要是因为其积沙量较多, 防沙效益明显(谢胜波等, 2014)。笔者的模拟结果表明, 栅栏沙障对风速的减弱作用更强, 影响范围也大于PE阻沙网, 而且栅栏沙障拦沙能力强于PE阻沙网, PE阻沙网输沙率较高, 障前几乎没有积沙, 障后沙粒的堆积高度较低, 更容易输送至路基处, 因此在布设防护措施时可适当增加栅栏沙障的数量, 提高防沙效益。李凯崇等(2015)对几种不同的挡沙墙进行了风洞试验研究, 其结果也表示轨枕式挡沙墙对风能的削弱能力远大于高立式PE网沙障, 并且PE网沙障输沙率大于其他高立式沙障。马学喜等(2016)发现不同结构的HDPE (High density PE)阻沙网都能阻拦大量流沙, 积沙主要分布于阻沙网的下风侧, 而上风侧积沙量相对较少。上述研究结果与本文模拟的结果近似, 说明模拟的结果基本可以反映路基及沙障的流场和积沙分布的实际情况。

本研究中, PE阻沙网后积沙范围在2H内, 并且阻沙网前几乎没有积沙, 这与郝晓杰等(2012)和张克存等(2016)的研究略有出入, 特别是障后积沙范围差别较大。郝晓杰的研究指出PE阻沙网后积沙范围在10H左右, 张克存等人(2016)利用Leica C10三维激光扫描仪精确测量了北麓河段两侧典型防沙措施及防护体系积沙形态特征和风沙堆积量, 其结果表明尼龙网阻沙栅栏障后积沙范围可达5H左右但障前并无积沙。本研究模拟的风速为固定值, 但实际情况下风速会不断变化, 风速较小的时候沙粒会堆积在靠近阻沙网的位置, 大风会把靠近阻沙网的沙粒继续向后输送并沉积, 这会导致PE阻沙网后积沙范围越来越大, 最终达到阈值。同时, 由于研究所选的材料及疏透度不同, 也可能会造成研究结果的不一致。

尽管数值模拟结果显示笔者设计的输沙铁路路基及沙障的组合形式明显减少了路基周围积沙量, 然而目前并未在青藏铁路实际采用, 尚无法实际验证其疏导效果, 需要在今后的研究中继续探索其优缺点, 在实践中进一步检验和完善。

4 结论

3种路基在无防护措施的情况下, 均有不同程度的积沙, 普通路基积沙情况最严重, 迎风坡积沙量明显大于背风坡; 通风路基边坡积沙较少, 部分沙粒沉积在通风管, 破坏通风管的保温功能; 输沙路基边坡几乎没有积沙, 并且输沙管内的积沙可以通过风力作用输移到积沙池内集中处理, 节省大量人力物力, 在实际应用中可以更好地避免风沙危害带来的损失。

在布设有沙障的情况下, 3种路基边坡的积沙量都明显减少, 沙障在一定程度上阻挡了风沙流, 减少风沙对路基的危害。沙障与路基之间存在一个低速区, 该区域速度两端低, 中间高, 根据这个规律, 可以适当调整沙障布设的位置, 使二者之间的过渡区域风速无法达到起沙风速, 保证沙障背风侧积沙不被再次起动, 提高沙障的防沙效益, 减轻风沙对路基的侵蚀。

模拟结果表明, 输沙路基相对于普通路基和通风路基更适合布设在青藏铁路风沙路段, 其独特的设计可有效地增加来流风速, 防止风沙上道, 并且有专门的积沙池收集堆积的沙粒, 但沙粒不能完全沉积在积沙池内, 会有少量沙粒沿积沙池攀升, 堆积在积沙池后。输沙路基的设计为青藏铁路防沙治沙提供了一种新的思路, 从研究结果来看, 此设计行之有效, 可在实践中进行尝试, 为防治沙害做出贡献。