呂帥元 田榮燕** 萬(wàn)常興 馬仁躍

1.西藏大學(xué)工學(xué)院，西藏 拉薩 850000；2.西藏自治區(qū)公路局青藏公路分局安多公路養(yǎng)護(hù)段，西藏 那曲 852000

風(fēng)沙流是沙質(zhì)地表上的沙塵顆粒受氣流影響產(chǎn)生地表蠕移、地表遷移、懸移并發(fā)生堆積的氣固兩相流現(xiàn)象［1］，其90%以上的沙粒徑在0.074～0.25mm 之間，主要發(fā)生在干旱或半干旱地區(qū)［2］。近年來，學(xué)者們從不同角度研究了風(fēng)沙流特征及其沙粒運(yùn)動(dòng)模式，如Long Shi［3］等人用數(shù)值模擬的方法研究了風(fēng)沙活動(dòng)特征、風(fēng)沙流粒度分布和沙粒濃度，并提出了整治風(fēng)沙危害的相應(yīng)措施；劉芳［4，5］運(yùn)用二維歐拉顆粒兩相流模型及一種新的宏觀顆粒模型（MPM）研究沙粒在氣流中的運(yùn)動(dòng)，得到近地面風(fēng)沙流的沙相速度分布，實(shí)現(xiàn)了風(fēng)沙兩相流的精確模擬；Lorenzo Raffaele［6］研究了不同風(fēng)向角、有效剪切速度和風(fēng)流特征等條件下的沙粒沉積模式，得出了以沙通量邊界條件和沙床初始狀態(tài)為變量的風(fēng)沙變化趨勢(shì)以及沙粒運(yùn)輸和堆積形態(tài)的動(dòng)力學(xué)演變；叢順［7］等人通過風(fēng)沙風(fēng)洞試驗(yàn)，模擬了實(shí)際環(huán)境下的流場(chǎng)分布特征，得出了流場(chǎng)中沙濃度梯度分布規(guī)律。

風(fēng)沙現(xiàn)象的出現(xiàn)不僅會(huì)對(duì)周圍生態(tài)造成破壞，也會(huì)對(duì)處于風(fēng)沙環(huán)境中的公路鐵路等基礎(chǔ)設(shè)施產(chǎn)生影響，沙粒沉積于路面，對(duì)行車安全及鐵路安全造成威脅。為此，眾多學(xué)者在了解風(fēng)沙運(yùn)動(dòng)模式的基礎(chǔ)上，針對(duì)路面積沙現(xiàn)象進(jìn)行了研究，如張帥［8］等人通過對(duì)公路護(hù)欄進(jìn)行挖孔改造和風(fēng)洞試驗(yàn)，分析了不同孔洞防護(hù)欄對(duì)路面風(fēng)速流場(chǎng)的影響，并計(jì)算了絕對(duì)含沙量，對(duì)未來不同行車環(huán)境的公路護(hù)欄改造提供了寶貴參考；Cui Wang［9］通過風(fēng)洞試驗(yàn)和現(xiàn)場(chǎng)試驗(yàn)?zāi)M了4種公路護(hù)欄附近的風(fēng)流場(chǎng)以及沙粒沉積特性，孫興林［10］等人通過Fluent 軟件研究了青藏鐵路原有路基與新型輸沙路基的流暢分布情況以及積沙特征，為青藏鐵路防沙治沙提供參考。

本文針對(duì)G109 國(guó)道部分路段存在的路面積沙現(xiàn)象，利用軟件COMSOL Multiphysics 6.0 對(duì)風(fēng)沙流運(yùn)動(dòng)進(jìn)行仿真模擬并結(jié)合計(jì)算流體動(dòng)力學(xué)理論得出路面風(fēng)速流場(chǎng)的分布狀況，對(duì)沙粒進(jìn)行運(yùn)動(dòng)追蹤得出其運(yùn)動(dòng)軌跡，并分析沙粒沉積特性，為處于風(fēng)沙環(huán)境中的公路護(hù)欄選型提供參考。

1 積沙路段環(huán)境數(shù)據(jù)分析

積沙路段位于G109北京—拉薩公路K3462～K3463之間，海拔4500m 以上，與正北向夾角42.08°，植被稀疏，土地沙化嚴(yán)重，研究區(qū)位置及現(xiàn)場(chǎng)積沙狀況如圖1所示。對(duì)該路段環(huán)境數(shù)據(jù)進(jìn)行監(jiān)測(cè)，包括環(huán)境溫度，風(fēng)速風(fēng)向等。通過對(duì)監(jiān)測(cè)數(shù)據(jù)的整理，提取出2021年5 月15 日至2022 年5 月15 日風(fēng)速有效數(shù)據(jù)103563組，選取每日10分鐘平均風(fēng)速的最大值及風(fēng)向作為統(tǒng)計(jì)量，結(jié)果如圖2所示。

圖1 研究區(qū)位置及積沙現(xiàn)狀圖

圖2 風(fēng)速分布圖

統(tǒng)計(jì)結(jié)果顯示，一年內(nèi)風(fēng)速平均值為7.8m/s，風(fēng)速處于5～10m/s 的占比79.23%，風(fēng)向以SSW 和SW 為主。結(jié)合積沙路段走向，風(fēng)沙流與公路夾角近乎垂直，可按垂直風(fēng)向進(jìn)行分析。

2 模擬方法

模擬分析需進(jìn)行如下假設(shè)：

（1）不考慮風(fēng)沙流對(duì)路段縱向的影響；（2）模型忽略立柱及連接細(xì)節(jié)，僅做截面分析；（3）風(fēng)速穩(wěn)定，不考慮濕度對(duì)沙流的影響；（4）沙粒為均勻大小的隨機(jī)分布顆粒。

2.1 幾何模型

幾何模型包括兩部分：路基和護(hù)欄。本文中選取積沙路段現(xiàn)有三波形梁鋼護(hù)欄和其他兩種常見護(hù)欄，即混凝土護(hù)欄和纜索護(hù)欄。

路基結(jié)構(gòu)模型按照積沙路段實(shí)際尺寸建模，相對(duì)高程3.2m，路面寬8m，迎風(fēng)側(cè)邊坡長(zhǎng)5m，背風(fēng)側(cè)邊坡長(zhǎng)4m。根據(jù)實(shí)地測(cè)量得知：三波形梁鋼護(hù)欄面板寬50cm，波峰高8.5cm，中心線高于路面70cm?；炷磷o(hù)欄與纜索護(hù)欄根據(jù)《公路交通安全設(shè)施設(shè)計(jì)細(xì)則》（JTG/T D81-2017）［11］進(jìn)行建模，均采用二（B）級(jí)防護(hù)等級(jí)，纜索護(hù)欄鋼索直徑均為1.8cm，四根鋼索間距13cm，最下一根高度43cm，；混凝土護(hù)欄高于路面81cm，頂面寬17.2cm，底面寬30cm。風(fēng)洞長(zhǎng)55m，高20m，模型及風(fēng)洞尺寸如圖3所示。

圖3 模型及風(fēng)洞尺寸平面圖

2.2 參數(shù)與邊界條件

選用標(biāo)準(zhǔn)k-ε瑞流模型，最大迭代次數(shù)100次，沙粒體積分?jǐn)?shù)為0.01。邊界均為無滑移流體壁面，左側(cè)為流體入口，右側(cè)邊界為流體出口，上邊界為開放邊界，下邊界及路基表面為沙粒沉積邊界，啟用粒子累積計(jì)數(shù)器，沙粒沉積邊界Y 向變形與粒子累積計(jì)數(shù)器相關(guān)聯(lián)來計(jì)算積沙厚度，計(jì)算時(shí)間30s 以反映沙粒運(yùn)動(dòng)軌跡及積沙的分布規(guī)律。綜合文獻(xiàn)研究結(jié)果［12，13］，模型主要計(jì)算參數(shù)如表1所示。

表1 計(jì)算參數(shù)

2.3 基本理論

2.3.1 流體運(yùn)動(dòng)方程。將風(fēng)沙流視為不可壓縮牛頓流體［14］，不考慮密度與運(yùn)動(dòng)黏度的變化，平面運(yùn)動(dòng)滿足方程：

（1）質(zhì)量守恒

其中：ρ為流體密度；▽為Nabla算子：

（2）動(dòng)量守恒

均質(zhì)、不可壓縮流體動(dòng)量方程為：

其中：μ為動(dòng)力粘度；P為運(yùn)動(dòng)壓力為體力。

（3）湍流

湍流動(dòng)能k守恒方程：

湍流動(dòng)能耗散率ε方程：

2.3.2 沙粒運(yùn)動(dòng)方程。不考慮沙粒的旋轉(zhuǎn)，單個(gè)沙粒在運(yùn)動(dòng)中受到重力及氣流曳力的作用［15］，其運(yùn)動(dòng)方程為：

與沙粒的尺寸、質(zhì)量和速度以及氣流的粘性、速度相關(guān)：

將重力及曳力公式代入式（5），則沙粒平面運(yùn)動(dòng)中X、Y 方向的運(yùn)動(dòng)方程為：

3 結(jié)果與分析

3.1 水平流場(chǎng)速度分布特征

風(fēng)速場(chǎng)是表示氣流運(yùn)動(dòng)的綜合指標(biāo)，沙的輸送量隨風(fēng)速的增大而增加，而速度較弱的區(qū)域通常會(huì)遇到沙子沉積問題［16］。事實(shí)上，自然風(fēng)場(chǎng)中的流動(dòng)大多為瞬態(tài)演變的，因此本文中采用瞬態(tài)求解器。流體流速初始值為0，為使流場(chǎng)入口邊界處初始值平滑過渡，且更接近自然風(fēng)場(chǎng)，本文中定義區(qū)間0～1的二階可導(dǎo)矩形波函數(shù)（wv），函數(shù)周期為5π，過渡區(qū)為3，X 軸與時(shí)間關(guān)聯(lián)，并賦予速度V=7.8［m/s］×wv(t［1/s］)，函數(shù)曲線如圖4。考慮流體慣性項(xiàng)，圖5反映了氣流遇到路基和護(hù)欄之后的風(fēng)速場(chǎng)演變過程。

圖4 矩形波函數(shù)曲線

圖5 水平流場(chǎng)速度分布

路基模型相同，迎風(fēng)側(cè)形成的初始流場(chǎng)分布具有相似性，迎風(fēng)側(cè)邊坡對(duì)水平氣流具有阻礙和抬升作用，形成半圓形速度減弱區(qū)，坡腳風(fēng)速接近0。邊坡對(duì)氣流的抬升作用使得路面斜上方流量增大，路肩出現(xiàn)風(fēng)速加速區(qū)，護(hù)欄種類的差異造成加速區(qū)的形成范圍大小不一，混凝土護(hù)欄形成的加速區(qū)范圍最大，纜索護(hù)欄最小。受護(hù)欄遮蔽作用的影響，氣流在路面上方的變化出現(xiàn)較大差異：波形護(hù)欄將流經(jīng)的氣流分離成兩部分，下部流經(jīng)面板與路面間的空心區(qū)域，上部通過護(hù)欄上方，且速度均有增加，氣流在護(hù)欄順風(fēng)處發(fā)生二次分流，在護(hù)欄背風(fēng)處形成逆流，與護(hù)欄下方空心處的正流相結(jié)合，逆流被抵消，正流被減弱，形成錐形渦流區(qū)；混凝土護(hù)欄由于不透風(fēng)性，將左側(cè)來流全部阻擋，護(hù)欄頂部邊緣處風(fēng)速上升較大，護(hù)欄背風(fēng)側(cè)形成大范圍三角形弱風(fēng)區(qū)；纜索護(hù)欄具有良好的通透性，路面上方風(fēng)速接近初始風(fēng)速，但仍會(huì)因?yàn)殇撍鞯姆至髯饔贸霈F(xiàn)加速區(qū)。路基背風(fēng)側(cè)均出現(xiàn)不同程度的渦流區(qū)，風(fēng)速低于2m/s，在距坡腳20m 處，風(fēng)速逐漸上升。由于矩形波狀來流的速度變化，中間時(shí)刻模型周圍流場(chǎng)分布較凌亂，流體受重力、慣性及邊界壓強(qiáng)的作用，在路面上方均出現(xiàn)不同程度的回流。

3 種護(hù)欄在兩側(cè)路肩0.5m 處25s 時(shí)風(fēng)速隨高度變化關(guān)系如圖6 所示。由圖6(a)可以看出：波形護(hù)欄與纜索護(hù)欄在迎風(fēng)側(cè)路肩0.5m 高度內(nèi)風(fēng)速變化趨勢(shì)相同，均為先增大后減小的趨勢(shì)，波形梁鋼護(hù)欄的峰值速度達(dá)到10.93m/s，纜索護(hù)欄的峰值速度為11.1m/s，由于受到纜索的影響，風(fēng)速出現(xiàn)鋸齒狀波動(dòng)；1m 高度以上，三種護(hù)欄附近的風(fēng)速均穩(wěn)定減小，10m 高度處風(fēng)速分別接近8.64m/s(波)，8.85m/s(混)和8.48m/s(纜)，風(fēng)速差值不超過0.5m/s，且在逐漸縮小，由此可知，在一定高度處，風(fēng)速流場(chǎng)的分布將不再受到護(hù)欄、路基的影響。圖6(b)為背風(fēng)側(cè)路肩的風(fēng)速曲線，波形梁鋼護(hù)欄與纜索護(hù)欄的風(fēng)速變化趨勢(shì)基本相同，但2m高度為一個(gè)分界點(diǎn)，2m以上的高度，風(fēng)速先增大后減小，波形梁鋼護(hù)欄在10m 高度處風(fēng)速接近9.17m/s，纜索護(hù)欄10m 高度處風(fēng)速接近8.64m/s；2m 以下的高度，風(fēng)速也是先增大后減小，但纜索護(hù)欄的峰值風(fēng)速明顯更高。

圖6 路肩風(fēng)速隨高度變化

3.2 沙粒運(yùn)動(dòng)軌跡

沙粒主要存在于0.5m高度內(nèi)的氣流中，故本文在氣流入口處設(shè)置0.5m高沙流，采用COMSOL流體流動(dòng)粒子追蹤模塊，下邊界采用變形網(wǎng)格，與粒子累積計(jì)數(shù)器建立聯(lián)系，沙粒運(yùn)動(dòng)受湍流場(chǎng)影響并考慮重力，運(yùn)動(dòng)軌跡如圖7所示。

圖7 沙粒運(yùn)動(dòng)軌跡

模型下邊界（包括路基邊界）為沙粒沉積邊界，不同顏色代表沙粒運(yùn)動(dòng)的速度，紅色最大，藍(lán)色最小。從圖7 可以看出，風(fēng)沙流在經(jīng)過路基護(hù)欄時(shí)運(yùn)動(dòng)軌跡并不相同：沙粒經(jīng)過混凝土護(hù)欄時(shí)會(huì)隨抬升的氣流出現(xiàn)較高的上揚(yáng)，高度可達(dá)10m，且運(yùn)動(dòng)距離較遠(yuǎn)，由于重力及回流作用，沙粒會(huì)出現(xiàn)反方向運(yùn)動(dòng)，直至沉積于背風(fēng)側(cè)坡腳；沙流經(jīng)過波形護(hù)欄時(shí)，面板上方與下方均有沙粒通過，根據(jù)對(duì)流場(chǎng)的分析可知，面板下邊緣處風(fēng)向與路面有一定夾角，且風(fēng)速較大，沙粒受此處氣流影響，傾斜運(yùn)動(dòng)至路面后半段，又因?yàn)橹亓帮L(fēng)速減弱的緣故出現(xiàn)沉積；沙流經(jīng)過纜索護(hù)欄時(shí)，其良好的通透性對(duì)沙粒運(yùn)動(dòng)的影響較小，沙粒在路面上方做近似拋物線運(yùn)動(dòng)，可完全躍過路面，在路基背風(fēng)側(cè)受回流與重力作用在坡腳出現(xiàn)沉積，但根據(jù)沙粒在路面上方躍起的高度與初始風(fēng)速可推斷，若風(fēng)速下降至起沙閾值，沙粒仍有很大可能沉積于路面。

3.3 積沙分布特征

對(duì)于波形護(hù)欄，沙粒首先沉積在迎風(fēng)側(cè)邊坡，隨著風(fēng)沙流的持續(xù)流動(dòng)，很容易在受到護(hù)欄干擾的路面上方沉積，路面積沙主要集中在右半段，與實(shí)際觀測(cè)到的情況比較相符。風(fēng)沙流運(yùn)動(dòng)25s 后，少部分沙粒落至背風(fēng)側(cè)邊坡及坡腳，厚度為迎風(fēng)側(cè)坡腳的37.2%.沙粒在迎風(fēng)側(cè)邊坡沉積的過程中，若遇到短暫且流速過大的氣流影響，會(huì)將沉積的沙粒再次卷上路表，進(jìn)一步加大路面積沙厚度。對(duì)于混凝土護(hù)欄，積沙主要分布于迎風(fēng)側(cè)邊坡、護(hù)欄根部和背風(fēng)側(cè)坡腳，路面也有少量積沙，造成這種分布的主要原因是沙粒運(yùn)動(dòng)的過程中受到邊坡的阻擋以及重力作用，分布范圍反而會(huì)受到壓縮，在這種情況下，僅有部分沙粒受氣流的抬升作用和沙粒間的碰撞作用躍過護(hù)欄，其余部分沙粒會(huì)在迎風(fēng)側(cè)坡腳與護(hù)欄根部沉積。對(duì)于纜索護(hù)欄，由于風(fēng)動(dòng)量基本保持不變，沙通量分布合理，因此沙粒順利通過路基。積沙主要分布于兩側(cè)坡腳，從圖8（c）可以看出，坡腳積沙寬度大致相同，但迎風(fēng)側(cè)坡腳積沙厚度比背風(fēng)側(cè)坡腳積沙厚度多1 倍。

圖8 積沙分布圖

以圖8 來看，三種護(hù)欄模型在背風(fēng)側(cè)坡腳處均有不同程度的積沙，這種情況可以理解為：飽和沙流通過路基的過程中，部分沙粒沉積，但風(fēng)速在水平方向上為瞬態(tài)演變，因此沙流可能變?yōu)椴伙柡蜕沉?，直到背風(fēng)側(cè)，風(fēng)速降低，與湍流相互影響，氣流變?yōu)轱柡蜕沉骰蜻^飽和沙流，沙粒再次沉積。因此，背風(fēng)側(cè)坡腳沙粒沉積在很大程度上與路基有關(guān)。

4 結(jié)論

（1）氣流在路基護(hù)欄模型迎風(fēng)側(cè)形成的流場(chǎng)分布具有相似性，隨著流場(chǎng)的發(fā)展，波形護(hù)欄面板將來流分為兩部分，面板邊緣處風(fēng)速增大；混凝土護(hù)欄背風(fēng)側(cè)形成較大三角形弱風(fēng)區(qū)；纜索護(hù)欄周圍流場(chǎng)分布簡(jiǎn)單，風(fēng)速接近初始風(fēng)速；中間時(shí)刻，來流速度為0m/s時(shí)，流體受重力、慣性及邊界壓強(qiáng)的作用，模型上方出現(xiàn)不同程度的回流；隨著高度的增加，三種護(hù)欄類型對(duì)流場(chǎng)的分布影響逐漸減小，由此可知，在一定高度處，風(fēng)速流場(chǎng)的分布將不再受到護(hù)欄、路基的影響。

（2）沙粒經(jīng)過混凝土護(hù)欄時(shí)會(huì)被氣流抬升，纜索護(hù)欄附近沙粒做類拋物線運(yùn)動(dòng)，一部分沙粒在波形梁鋼護(hù)欄路面背風(fēng)側(cè)貼地滑行；由于重力及來流風(fēng)速的變化，沙粒會(huì)隨著回流作反方向運(yùn)動(dòng)，容易造成沉積。

（3）三種模型在路基迎風(fēng)側(cè)與背風(fēng)側(cè)邊坡均有不同程度的積沙，只有纜索護(hù)欄在路面上方積沙較少，波形護(hù)欄路面積沙主要分布在右半段；沙流運(yùn)動(dòng)25s后，三種護(hù)欄模型背風(fēng)側(cè)邊坡積沙厚度約為迎風(fēng)側(cè)邊坡的一半。