張?zhí)旄?，周晅毅，顧 明

(同濟大學土木工程防災(zāi)國家重點實驗室，上海 200092)

在我國北部的多雪地區(qū)，橋面積雪會導(dǎo)致車輛難以通行，致使道路受阻；由風引起的橋面不均勻分布雪荷載可能會加劇該交通問題。為減小橋面風致積雪對交通的不利影響，需要對橋面雪飄移的機理進行探究，合理地預(yù)測橋面積雪的重分布形式，進而采取相應(yīng)的措施。由于橋梁附近的流場受護欄的影響較大，在探究橋面雪飄移時需要考慮護欄的影響。

雪飄移問題的研究方法主要有現(xiàn)場實測[1]、風洞試驗[2]和數(shù)值模擬[3]。近年來，隨著計算機技術(shù)的發(fā)展，基于計算流體動力學(CFD)的數(shù)值模擬方法逐漸興起，成為工程領(lǐng)域中研究雪飄移問題的有效手段[4]，并廣泛用于地面或建筑屋面雪飄移的模擬[5-8]。TOM INAGA 等[7]對建筑周邊雪飄移的CFD模型進行了概述，并提出了一個模擬雪飄移的新模型。ZHOU 等對平屋面[8]、雙坡屋面[9]上的雪飄移進行了模擬，預(yù)測了屋面上的風致積雪重分布。ZHU 等[10]對雪飄移的準定常模擬方法進行了改進，并將其用于高低屋面雪飄移的研究。周晅毅等[11]采用拉格朗日方法，對立方體周邊的雪飄移進行了模擬，并考慮了顆粒運動后的質(zhì)量濃度對雪層表面摩擦速度的修正。MA 等[12]結(jié)合現(xiàn)場實測和數(shù)值模擬，對路堤周圍風致積雪堆積進行了預(yù)測。

然而，當前針對橋面雪飄移問題的研究依然較少。DUAN 等[13]采用數(shù)值模擬方法對橋面風致積雪重分布進行了預(yù)測，但其主要關(guān)注橋面的動力特性，沒有對不同護欄下橋面雪飄移的特點進行探討。

本文基于CFD方法，探究不同護欄透風率下橋面雪飄移的特點。首先通過一個平屋面案例，對本文數(shù)值模擬方法的正確性進行了驗證。隨后采用該方法進行橋面雪飄移的數(shù)值模擬。橋面上共設(shè)置4 道防撞欄，在模擬的過程中改變護欄的透風率(0%～100%)，對比不同護欄透風率下橋面附近流場形式、摩擦速度、橋面雪濃度和橋面風致積雪重分布的異同，基于模擬結(jié)果對實際工程中橋梁護欄的設(shè)計給出建議。

1 數(shù)值模擬方法

1.1 數(shù)值模擬方案

風致積雪運動屬于兩相流運動，本文采用定常的歐拉-歐拉方法進行模擬。空氣相采用Reynolds時間平均方法，湍流模型采用realizablek-ε 模型；雪相采用歐拉方法，通過附加一個雪相輸運方程來模擬。模擬時假設(shè)空氣相和雪相之間為單向耦合的，即雪在風的作用下發(fā)生飄移，而雪的飄移對風場不產(chǎn)生影響[14]。

1.2 雪相控制方程

根據(jù)雪顆粒運動高度的不同，可以將雪的運動形式劃分為：蠕移(小于0.01 m)、躍移(0.01 m～0.1 m)和懸移(0.1 m～100m)，在模擬時通常將蠕移包含在躍移中[15]。在懸移運動時，雪顆粒會因重力作用產(chǎn)生向下的沉降；而在躍移運動時，由于雪顆粒會在雪表面發(fā)生回彈，重力沉降的作用會被減弱。雪相的控制方程為[16]：

式中： φ為雪濃度；νt為湍流運動粘度；Sct為湍流施密特數(shù)；wf為沉降速度，在躍移層(wf,sal)和懸移層(wf,sus)取值不同。等式左側(cè)分別為瞬態(tài)項和對流項，等式右側(cè)分別為擴散項和考慮重力沉降的附加對流項。由于本文進行定常的數(shù)值模擬，因此不考慮雪隨時間的瞬態(tài)變化(即瞬態(tài)項為0)，只考慮雪的對流和擴散。

1.3 雪的侵蝕和沉積

雪表面的摩擦速度是衡量雪侵蝕和沉積的關(guān)鍵因素，其表示雪表面所受剪切力的大小，定義為：

式中：τ0為風作用下雪表面的剪切力；ρa為空氣的密度。

雪的侵蝕和沉積通量可以采用NAAIM 等[17]提出的經(jīng)驗公式進行計算：

式中：Aero=0.0007為經(jīng)驗系數(shù)；u*t為閾值摩擦速度；下標ero表示侵蝕；acc表示沉積?？梢钥闯?，當雪表面的摩擦速度大于閾值摩擦速度時，雪表面會發(fā)生侵蝕，反之則產(chǎn)生沉積。

基于雪侵蝕和沉積通量的計算結(jié)果，可以得到單位時間積雪高度的變化量：

式中，ρb為雪的堆積密度。

2 數(shù)值模型

2.1 研究對象

本文采用一個典型的分體箱梁斷面作為研究對象[18]，橋面寬度(L)為49.0 m，斷面高度(H)為4.0m，橋面初始積雪厚度為0.3m，橋面兩側(cè)和中部共設(shè)置4道防撞欄(以下簡稱“護欄”)，護欄的高度為1.6m(高出雪面部分)，根據(jù)護欄位置將橋面劃分為5 個區(qū)域，如圖1所示。護欄的透風率從100%(無護欄)到0%(不透風護欄)變化，變化間隔為12.5%，根據(jù)護欄透風率的大小可以將其分為：無護欄、高透風率護欄和低透風率護欄三種類型，見表1。

表1 工況描述Table1 Description of cases

圖1 橋梁斷面/cm Fig.1 Bridge section profile

2.2 數(shù)值模擬計算域和網(wǎng)格劃分

由于橋梁軸向長度大，斷面形狀沿軸向基本保持不變，因此對橋梁斷面建立二維模型。利用ANSYS/FLUENT軟件對橋面附近的風場進行模擬；將雪相的輸運方程作為額外的方程對軟件進行二次開發(fā)，從而實現(xiàn)雪飄移的模擬。圖2為橋面數(shù)值模擬的計算域。網(wǎng)格方案采用漸變的結(jié)構(gòu)化網(wǎng)格，網(wǎng)格總數(shù)為7.0萬，在橋面附近對網(wǎng)格進行了加密處理。圖3為各種護欄透風率工況下橋面附近的網(wǎng)格，橋面首層網(wǎng)格的高度為0.05m，能夠?qū)S移層和懸移層進行區(qū)分。經(jīng)檢驗，在該網(wǎng)格方案下，計算結(jié)果具有較好的網(wǎng)格無關(guān)性。

圖2 數(shù)值模擬計算域Fig.2 Computational domain of numerical simulation

圖3 橋面附近網(wǎng)格Fig.3 Grid around bridge deck

2.3 數(shù)值模擬參數(shù)

在數(shù)值模擬的過程中，需要對模擬參數(shù)進行合理的選取[9]。我國北方天氣干燥，雪顆粒易在風的作用下發(fā)生飄移，且在橋面上發(fā)生飄移的雪大多是新降雪，于是本文考慮干燥、新降雪的情況，雪的堆積密度取70 kg/m3[19-21]，閾值摩擦速度取0.07 m/s[22-23]。湍流施密特數(shù)取1.0[15]。懸移層受重力沉降的影響較大，其沉降速度取0.5m/s[15]，而躍移層受重力沉降的影響相對較小，其沉降速度取0.15 m/s。初始積雪深度取30 cm[24]。數(shù)值模擬的參數(shù)匯總見表2。

表2 雪飄移的數(shù)值模擬參數(shù)匯總Table2 Parameters in numerical simulation of snow drifting

2.4 邊界條件設(shè)置

來流風速在橋梁高度方向上的變化不顯著，因此在速度入口采用均勻的風剖面。根據(jù)氣象資料，對我國北部某地區(qū)冬季風速數(shù)據(jù)進行廣義極值分布(GEV)擬合，如圖4所示，得到具有95%保證率的估計風速：8.5m/s。于是，在本文的數(shù)值模擬中，取入口處的風速為8.5m/s，入口處的湍動能k和湍流耗散率ε 用下式計算[25]：

圖4 我國北部某地區(qū)冬季風速數(shù)據(jù)(1980-2014)Fig.4 W ind velocity during w inter in representative northern region of China (1980-2014)

式中：I為湍流強度，本文取0.12；Cμ=0.09為常數(shù)；μ為湍流動力粘度；μt/μ為湍流粘度比。

速度出口邊界設(shè)置為自由出流邊界，上下邊界設(shè)置為對稱邊界，雪表面采用標準壁面函數(shù)。湍流模型采用realizablek-ε 模型，離散格式為二階迎風格式。收斂標準為無量綱殘差降至10-7以下，且收斂時橋面附近監(jiān)測點的風速保持不變。

3 試驗驗證

3.1 驗證對象

由于平屋面和橋梁斷面在外形上具有一定的相似性，二者表面雪輸運的條件也相似，因此選擇平屋面作為本文的驗證對象。以文獻[26]中的平屋面為驗證對象，該屋面原型高(HR)6.0m，長(LR)36.0m。

3.2 平屋面雪飄移風洞試驗介紹

ZHOU 等[26]在同濟大學1號風洞進行了平屋面積雪重分布的風洞試驗。風洞試驗的縮尺比為1∶25，縮尺后的試驗?zāi)Ｐ烷L144 cm，高24 cm，如圖5所示。該試驗采用硅砂模擬雪顆粒，屋面模型上的初始硅砂厚度為4.0 cm，對應(yīng)原型的初始積雪深度為1.0m。風洞試驗的阻塞比滿足小于5%的要求。限于篇幅，風洞試驗的詳細介紹見文獻[26]。

圖5 平屋面風洞試驗?zāi)Ｐ?cm Fig.5 Flat roof model in w ind tunnel test

3.3 平屋面雪飄移數(shù)值模擬介紹

以文獻[26]中的平屋面原型為模擬對象，進行平屋面雪飄移的數(shù)值模擬。數(shù)值模擬的計算域和網(wǎng)格劃分見圖6。為使該數(shù)值模擬條件與試驗的條件一致，將初始積雪深度設(shè)置為1.0 m，其余數(shù)值模擬參數(shù)的設(shè)置與橋面雪飄移模擬中的一致。

圖6 平屋面雪飄移數(shù)值模擬計算域和網(wǎng)格劃分Fig.6 Computational domain and grid in numerical simulation of snow drifting around flat roof

3.4 驗證結(jié)果

根據(jù)1.3節(jié)所述雪侵蝕沉積的計算方法，可以計算出平屋面上各個位置處積雪深度的變化，進而得到平屋面風致積雪的重分布。圖7為風洞試驗[26]和數(shù)值模擬積雪重分布的結(jié)果，其中，橫坐標x/L表示屋面上點的相對位置，縱坐標為無量綱的積雪深度，為屋面初始積雪深度(1.0 m)。從圖中可以看出，數(shù)值模擬結(jié)果和風洞試驗結(jié)果總體吻合較好，僅在部分區(qū)域存在差異：在x/L= 0.3附近，數(shù)值模擬結(jié)果比風洞試驗結(jié)果大8.6%；在x/L= 0.9附近，數(shù)值模擬結(jié)果比風洞試驗結(jié)果小6.5%。

圖7 平屋面積雪重分布Fig.7 Snow redistribution on flat roof

4 數(shù)值模擬結(jié)果

4.1 流場特性

圖8為橋面各區(qū)域的流線圖，僅給出護欄透風率為100%、50%和0%的三種工況，分別代表了無護欄、高透風率護欄和低透風率護欄的三種類型?？偟膩砜?，在各工況下，來流均在橋面的迎風前緣發(fā)生了分離，使橋面的迎風端產(chǎn)生了一個漩渦；在橋面的尾緣也產(chǎn)生了一個尾流渦。三種工況的差異主要體現(xiàn)在橋面中部(尤其是護欄附近)的區(qū)域。當護欄的透風率為100%和50%時，流場的形式大致相似；而當護欄的透風率為0%時，在護欄1的前端有另一個漩渦產(chǎn)生，此外，各護欄的背風側(cè)也都產(chǎn)生了漩渦，其中區(qū)域2(即護欄1的背風側(cè))中的漩渦尺寸最大。

圖8 橋面各區(qū)域流線Fig.8 Stream lines in each region on bridge deck

4.2 橋面(雪表面)摩擦速度

橋面摩擦速度是判斷雪侵蝕和沉積的重要參數(shù)，其與橋面附近的流場密切相關(guān)。圖9繪出了橋面摩擦速度的分布，相同類型的工況繪制在同一張圖上。總的來看，當護欄的透風率大于50%時(圖9(a)、圖9(b))，橋面摩擦速度相對較大，基本均大于閾值摩擦速度(0.07 m/s)；而當護欄的透風率小于50%時(圖9(c))，橋面摩擦速度相對較小，在某些區(qū)域小于閾值摩擦速度。

在高透風率護欄情況下，橋面摩擦速度隨護欄透風率的減小而減小，這體現(xiàn)了護欄對風的阻礙作用。而在低透風率護欄情況下，由于其流場形式復(fù)雜，摩擦速度隨透風率的變化趨勢也與高透風率護欄情況不同。以區(qū)域2為例，結(jié)合低透風率護欄情況下的流線圖(圖10)分析可知，當護欄透風率較小(但不為0%)時，護欄背風側(cè)會因狹管效應(yīng)形成一個氣流加速區(qū)(圖10)，導(dǎo)致護欄背風側(cè)的摩擦速度顯著增大(圖9(c))。此外，區(qū)域2的后部會出現(xiàn)一個漩渦，該漩渦的尺寸隨護欄透風率的減小而增大(圖10)。對比圖9(c)可以看出，漩渦中部的摩擦速度較大，邊緣的摩擦速度較小，且摩擦速度的極大值隨漩渦尺寸的增大而增大。

圖9 橋面摩擦速度Fig.9 Friction velocity on bridge deck

圖10 橋面區(qū)域 2 流線圖 (低透風率護欄情況)Fig.10 Streamlines in region No. 2 on bridge deck (cases with low barrier porosities)

為進一步分析橋面摩擦速度隨護欄透風率的變化，將橋面摩擦速度分別按各區(qū)域和全橋面進行平均，得到平均摩擦速度，如圖11所示。可以看出，橋面各區(qū)域平均摩擦速度都大于閾值摩擦速度(0.07m/s)，這意味著橋面積雪在強風的作用下主要發(fā)生侵蝕。在高透風率護欄情況下，橋面大多數(shù)區(qū)域的平均摩擦速度隨護欄透風率的減小而顯著減小，僅區(qū)域1(迎風側(cè)無護欄)的變化較為平緩。在低透風率護欄情況下，盡管護欄周圍存在加速區(qū)和漩渦區(qū)，但其影響范圍相對較小。從圖9(c)和圖11可以看出，低透風率時僅區(qū)域2的(平均)摩擦速度隨護欄透風率的減小而增大，其他區(qū)域摩擦速度隨透風率的變化不顯著。

圖11 平均摩擦速度Fig.11 Average friction velocity

4.3 橋面附近雪濃度

本文采用橋面首層網(wǎng)格雪濃度計算雪的沉積通量，首層網(wǎng)格的相對雪濃度見圖12，縱坐標φ / φmax表示用最大雪濃度無量綱化后的相對雪濃度。可以看出，橋面雪濃度和橋面摩擦速度的分布規(guī)律相似。這是因為當橋面摩擦速度較高時，雪顆粒被風挾卷到計算域中，經(jīng)對流和擴散后，該處雪濃度升高；而當橋面摩擦速度較低時，雪顆粒難以被風吹起，該處雪濃度也相應(yīng)較低。

圖12 橋面首層網(wǎng)格相對雪濃度Fig.12 Relative snow concentration at the first cell on bridge deck

4.4 橋面積雪重分布

圖13為2 d 后橋面積雪重分布的計算結(jié)果，縱坐標S/S0表示無量綱積雪深度，其中，S0= 0.3m為橋面初始積雪深度。在高透風率護欄(大于50%)情況下(圖13(a)、圖13(b))，橋面各區(qū)域僅發(fā)生雪的侵蝕，基本不發(fā)生雪的沉積，且橋面積雪的侵蝕程度隨護欄透風率的減小而減小，這體現(xiàn)了護欄對雪飄移的阻礙作用。同一工況下，區(qū)域1的侵蝕程度較小，而其他區(qū)域侵蝕程度相對較大，且各區(qū)域侵蝕的特征相似。在低透風率護欄(小于50%)情況下(圖13(c))，受護欄背風側(cè)加速區(qū)和漩渦區(qū)的影響，部分區(qū)域的積雪侵蝕程度隨護欄透風率的減小而增大，這與橋面摩擦速度的規(guī)律是一致的。低透風率時較為顯著的沉積主要發(fā)生在護欄1的背風側(cè)(區(qū)域2 漩渦邊緣)，該處摩擦速度小于閾值摩擦速度，因此雪顆粒不能被風挾卷進入計算域，只能產(chǎn)生沉積。隨著護欄透風率的減小，漩渦的尺寸增大(圖10)，沉積區(qū)域越來越靠近護欄(圖13(c))，最大積雪沉積量約為初始積雪深度(30 cm)的1.25倍。

圖13 橋面積雪重分布(2 d)Fig.13 Snow redistribution on bridge deck (2 days)

為了進一步分析橋面積雪隨護欄透風率的變化趨勢，將橋面積雪變化分別按各區(qū)域和全橋面進行了平均，得到平均積雪分布，如圖14所示。可以看出，在高透風率護欄情況下，橋面積雪的平均改變量隨護欄透風率的減小而減??；在低透風率護欄情況下，僅區(qū)域2的積雪平均改變量隨護欄透風率的減小而增大，其他區(qū)域的變化不顯著，這與平均摩擦速度(圖11)的規(guī)律是一致的。

圖14 平均積雪分布(2 d)Fig.14 Average snow redistribution (2 days)

在實際工程中，橋面上積雪的堆積會使道路受阻，為了減小風致積雪堆積對橋面交通的不利影響，建議在橋梁設(shè)計時采用高透風率(大于50%)的護欄。

5 結(jié)論

本文采用CFD方法，對橋面的風致積雪飄移進行了數(shù)值模擬，對比分析了不同護欄透風率下橋面風致積雪重分布形式，得出了以下幾個結(jié)論：

(1)以平屋面為案例，對比了屋面積雪重分布的CFD模擬結(jié)果和風洞試驗結(jié)果，兩種結(jié)果吻合較好，表明采用CFD模擬得到的結(jié)構(gòu)表面雪飄移結(jié)果可以滿足工程精度的要求。

(2)護欄透風率對橋面積雪分布形式有顯著的影響。當護欄透風率大于50%時，橋面上不會出現(xiàn)積雪的堆積，且橋面積雪侵蝕程度隨護欄透風率的減小而減??；同一工況下，區(qū)域1的侵蝕程度較小，而其他區(qū)域的侵蝕程度相對較大。當護欄透風率小于50%時，積雪會在橋面護欄附近產(chǎn)生沉積，且在迎風端護欄的背風側(cè)(區(qū)域2)沉積最嚴重，最大積雪沉積量約為初始積雪深度(30 cm)的1.25倍。

(3)基于數(shù)值模擬結(jié)果，建議在橋梁設(shè)計時采用高透風率(大于50%)的護欄，以減小風致積雪對橋面交通的不利影響。