復雜山區(qū)地形滑雪場區(qū)域風環(huán)境數值模擬

賀佳偉,趙亞哥白,張洪福,辛大波

(東北林業(yè)大學 土木工程學院, 黑龍江 哈爾濱 150040)

0 引言

滑雪場是提供高山滑雪、跳臺滑雪、越野滑雪等競技滑雪運動的體育場地,通常有賽道、運輸索道、服務場館等基礎設施。風場條件是運動員成績和滑雪場安全運營的主要影響因素之一,主要體現在2個方面。首先,運動員賽訓成績受到風速和風向的影響。其次,運輸索道的運行也受到風場環(huán)境的限制?；﹫鼋ㄔO選址大多數為山區(qū),滑雪場區(qū)域山丘起伏,溝壑縱橫,增加了滑雪場區(qū)域風場的復雜性,使得地處滑雪場區(qū)域的賽道、索道等基礎設施所受風的影響具有較強的不一致性。因此,開展復雜山區(qū)地形滑雪場區(qū)域風環(huán)境研究,精確地評價滑雪場區(qū)域的風環(huán)境,對于提高賽訓質量,保障安全運營具有重要意義。

開展復雜地形風場研究手段通常包括現場實測、風洞試驗和CFD(computational fluid dynamics, CFD)數值模擬?，F場實測是獲取風場信息的最有效手段[1-4],可以準確地分析風速譜、湍流強度、陣風因子、湍流積分尺度等脈動風特性。例如黃國慶等[5]通過在宣威方向橋位處附近建設測風塔,安裝數據采集系統(tǒng),實現了普立特橋位處風特性的實測研究。風洞試驗通過制作真實地形的縮尺模型,根據試驗目的合理布置監(jiān)測點來獲取監(jiān)測位置的風特性參數。例如YAMAGUCHI等[6]利用風洞試驗獲取平均風速和湍流度等,研究了日本沿海地區(qū)復雜地形下的氣流特性。胡峰強[7]以貴州北盤江大橋和湖北四渡河大橋周邊的地形為工程實例,結合風洞試驗結果以及附近氣象臺的數據資料,提出了復雜地形橋址地區(qū)設計基本風速梯度的修正方法。由于現場實測成本較高,需要長時間的數據采集,對儀器設備野外工作的性能有很強的要求[8],而風洞試驗受實驗室條件的限制,地形模型縮尺比較大時,會出現雷諾數效應,CFD數值模擬具有省時、經濟、高效的特點,因此,CFD數值模擬技術成為研究復雜地形風場特性的重要手段[9-11]。隨著CFD技術在復雜地形風場模擬的廣泛應用,在各方面的研究已取得很多進展。在復雜地形建模方面,胡朋等[12]和劉志文等[13]通過建立地形過渡段來消除來流風經過地形截斷邊緣出現的分離和繞流等現象,減小“人為峭壁”對數值模擬的結果影響。在模擬區(qū)域范圍大小及網格劃分方面,張希斌[14]研究了復雜地形范圍大小對數模擬的影響。肖儀清等[15]、周志勇等[16]和CHAO等[17]通過設計多種網格方案來確定最優(yōu)網格最優(yōu)劃分方式。在數值模湍流模型的選取方面,TANG等[18]采用開源CFD軟件OpenFOAM, 應用基于雷諾平均 Navier-Stokers(RANS)方程和修正的k-ε湍流模型,實現了高分辨率的風資源分布和準確的風速估算。HUANG等[19]基于平衡大氣邊界層模擬獲取入口邊界條件,用SSTk-ω模型對香港小蠔彎觀測站周圍的復雜地形進行多風向數值模擬,并與風洞試驗和觀測數據進行對比,驗證數值模擬的準確性。

綜上,近年來開展復雜地形風環(huán)境研究主要集中在工程建設選址、風能資源評估等方面,開展滑雪場區(qū)域風環(huán)境的研究較少。文中利用數值模擬方法針對亞布力滑雪場所在區(qū)域開展風環(huán)境研究,分析了滑雪場賽道區(qū)域、索道位置風場特點,并對賽訓和索道運營安全提出建議。

1 數值計算模型

1.1 工程概況

黑龍江亞布力滑雪場地處于山地林區(qū),是我國重要的滑雪比賽及訓練基地,雪場及設施建設具有典型的代表性。針對滑雪場區(qū)域地形的特點,從地形圖上選取滑雪場中心半徑為3 km范圍的地形。在此地形中包括1.2 km滑雪賽道和長800 m運輸索道,A～E為現場風速儀的安裝位置,滑雪場區(qū)域平面如圖1所示。雪場周圍具有數量眾多的高海拔山峰,風場環(huán)境具有高度的復雜性,滑雪場周邊高程如圖2所示,M1～M4為周邊山體,E為風速儀位置。

圖1 滑雪場區(qū)域平面圖Fig. 1 Plan of ski resort areas

圖2 滑雪場周邊高程圖Fig. 2 Elevation map around the ski resort

1.2 計算模型及網格劃分

針對圖1的地形范圍,高程數據采用全球數字高程模型(ASTER GDEM),分辨率為30 m?？紤]到地形數據直接建模會出現截斷高差,通過HUANG等[20]提出的地形過渡曲線建立過渡段來減小“人為峭壁”的影響,如圖3所示。過渡曲線如式(1):

圖3 地形過渡段Fig. 3 Terrain transition surfaces

(1)

式中,x、y為歸一化處理之后的過渡曲線的長度和高度,根據式(1)通過地形邊緣的高程就可以得到地形過渡段的長度。

為了確保來流風得到充分發(fā)展,根據FRANKE等[21]對數值模擬計算域的建議,入口距離地形中心的距離為5.5D,地形中心距出口距離為10.5D, 地形中心距兩對稱邊界的距離為5.5D,D是模擬地形的直徑6 km,計算域的高度Hd為7 km。計算域的尺寸為16D×11D×Hd, 地形模型阻塞率小于3%,滿足規(guī)范要求。

整個計算域分成3個部分,包括:外圍區(qū)域、空心圓柱過渡區(qū)域和內部圓柱地形區(qū)域,如圖4所示。

圖4 計算域網格劃分Fig. 4 Computational domain meshing

采用四面體和六面體網格混合的劃分方法,外圍區(qū)域由于關注來流和尾流發(fā)展的階段,選擇結構化六面體網格。過渡區(qū)域及地形區(qū)域,由于曲面不規(guī)則,采用非結構化的四面體網格。考慮到地形的復雜性,在近地面邊界層設置5層棱柱體網格,第1層網格高度為1 m,增長因子為1.1。Uf為入口10 m高度位置的風速(m/s),U為測點距離地面高度的風速(m/s),圖5為877萬、628萬、533萬的3種不同數量網格下,風速儀E位置的無量綱化風速剖面模擬結果。3種不同數量網格下風剖面幾乎完全重合,這表明文中計算結果與網格數量無關。為節(jié)約計算資源,在文中均選取533萬的網格模型進行滑雪場風環(huán)境的數值模擬研究。

圖5 網格無關性驗證Fig. 5 Independence analysis of grid resolution

1.3 邊界條件與計算參數設置

計算域內地形的粗糙度會對風速的垂直分布產生較大的影響,它在很大程度上決定了區(qū)域內的局部流動條件和近地面邊界層的發(fā)展[22]。Fluent在處理壁面函數時,是通過將粗糙度長度z0轉換為粗糙度物理高度ks進行設置[23],粗糙度長度和粗糙度物理高度之間存在的關系為:

(2)

式中:Cs為粗糙度常數,通常取0.5;根據地表粗糙度分類的劃分[24-25],有少量的樹或建筑的城鎮(zhèn),本研究中地形區(qū)域的地表粗糙度長度z0設為0.03 m。

計算域的上方和2個側面設置為對稱邊界條件,底面過渡區(qū)域及地形采用無滑移邊界條件,出口為壓力出口(pressure-out)。計算域入口通過用戶自定義函數(UDF)進行編譯,速度入口表達式為式(3),按我國規(guī)范的大氣邊界層B類地表風剖面規(guī)律(指數α為0.15)進行設置,湍動能k及耗散率ω的入口條件[20]如式(4)、式(5)所示:

(3)

(4)

(5)

式中:Ux表示地面高度z處的風速(m/s);Ur表示參考高度zr處的風速;湍流強度I取5%,湍流積分尺度l取500 m; c為常數0.033。

湍流模型選用對流動分離解析度較好的k-ωSST剪切應力模型[26],壓力與速度耦合處理選SIMPLEC算法。梯度插值方法為Green-Gauss Node Based,動量、湍動能、耗散率的控制離散格式采用二階迎風(second order upwind),壓力基非穩(wěn)態(tài)求解器,計算時間步長為0.1 s,總體殘差設為10-6,當速度迭代不再隨時間波動,認為計算收斂。

1.4 計算工況與觀測點設置

選取12個來流風向(見圖6),定義0°為正北風向,來流風的角度順時針旋轉為正,每30°為1個工況,工況1～12的來流風向為0°～330°?；﹫鲑惖绤^(qū)域地勢處于地形的凹段,選取風速儀E海拔高度處的平面風速,研究不同來流風向下地形對滑雪場賽道區(qū)域風環(huán)境的影響。

圖6 風向工況示意Fig. 6 Sketch of wind direction cases

為了得到索道處風環(huán)境特性,在索道纜車行進方向上站與下站之間等距布置9個觀測點,觀測點距地表高度約10 m,索道下站到上站地勢逐漸升高,上站與下站高程差約200 m,如圖7所示。 通過計算不同來流風下風速放大系數、風攻角、風偏角來獲取索道行進方向各觀測位置處的風環(huán)境特性。風速放大系數定義為山體某一高度風速的平均值與入口處同一高度處風速的比值,用以量化山體對來流風速的放大作用。

圖7 索道觀測點布置示意Fig. 7 Sketch of ropeway observation points layout

風速放大系數Cu、風攻角α、風偏角β的計算表達式見式(6)～式(8):

(6)

(7)

(8)

式中:U為索道觀測點位置的風速;Uf為入口10 m高度位置的風速;Uu、Uv、Uw、分別代表縱向、橫向、豎向的風速分量。

1.5 CFD數值結果驗證

基于前述的網格及相應的邊界條件,在山體前選取某一位置來繪制風速剖面,對來流風進行了自保持分析,從圖8中可以看出,山前與入口風剖面幾乎重合,說明風速具有較好的自保持性。

圖8 風剖面自保持性驗證Fig. 8 Self-preservation verification of wind speed profile

圖9給出了來流風經過地形過渡段后的流場情況,來流風在地形過渡段邊緣風速開始減小,經過地形過渡段后氣流有一定的緩慢抬升,當氣流到達地形中間區(qū)域時風速已基本保持來流風速的大小,表明地形過段曲線應用于山區(qū)風場建模具有較好的過渡性能,數值模擬結果具有可信度。

圖9 地形過渡區(qū)風速云圖Fig. 9 Wind speed contour of the terrain transition areas

為了驗證CFD數值模擬的準確性,根據圖1中各風速儀采集的數據,將CFD數值模擬結果與現場實測的風速及風向進行對比。定義2組風速比:

1)定義標準風速比K,K是測點位置的風速與參考點位置風速的比值。由于風速儀C的安裝位置海拔較高,風速受周圍地形地勢的影響較小,選取風速儀C的位置作為參考點,如式(9)、式(10)所示:

(9)

(10)

式中:Kexp、Knum為實測風速比和模擬風速比;Wi、Vi為風速儀i位置的實測和數值模擬的平均風速;WC、VC為風速儀C位置的實測和數值模擬的平均風速。

2)定義水平風速比KH,KH是某海拔高度處局部水平風速與模擬入口10 m高度處風速的比值,如式(11)所示:

(11)

式中:VH為海拔高度為H的局部水平風速;V10為遠方來流入口10 m高度處的平均風速。

觀測點C一個月(2019年12月)的現場實測風向統(tǒng)計結果如圖10所示,平均風向大多在90°～240°?，F場實測數據選取10 min平均風速數據,6個風向下的模擬風速比與實測風速比如圖11(a)所示。

圖10 風速儀C現場實測風向分布頻度Fig. 10 Frequency of wind direction distribution measured on site by anemometer C

從圖11(a)中可以看出,模擬風速比Knum與實測風速比Kexp的誤差基本在20%以內,表明模擬與實測對應方向的風速比具有良好的一致性。圖11(b)為數值模擬風向(PHInum)對應的現場實測方向(PHIexp),模擬風向與現場實測的偏差小于15°。

圖11 數值模擬結果與現場實測對比Fig. 11 Comparison between numerical simulation results and field measurements

來流風向分別為90°和180°,海拔高度為532 m的風速儀位置周邊區(qū)域的模擬風速比如圖12和圖13所示。入口風向為90°時來流風從開闊地帶向山谷地區(qū)靠近,風速逐漸減小。來流風向為180°時研究區(qū)內風速會受到高山影響,風速先減小后增加。數值結果符合該研究區(qū)地形風場的常規(guī)分布特征,表明數值模擬具有一定的準確性。

圖12 風速儀位置周邊區(qū)域KH云圖(來流風向90°,海拔高度532 m)Fig. 12 KH contour of the surrounding areas of the anemometer (incoming wind direction 90°, altitude 532 m)

圖13 風速儀位置周邊區(qū)域KH云圖(來流風向180°,海拔高度532 m)Fig. 13 KH contour of the surrounding areas of the anemometer (incoming wind direction 180°, altitude 532 m)

2 雪場區(qū)域數值模擬的應用

2.1 滑雪賽道風場分析

E點周圍的風速場分布如圖14所示,圖中黑色箭頭表示風的流動方向,云圖中地形范圍白色區(qū)域是山體在海拔高度為505 m處的截面。

圖14 E點周圍的KH云圖Fig. 14 KH contour of the surrounding areas of the point E

當來流風向為0°時,來流風在山體M1和M4之間出現加速,之后風速開始衰減,在到達賽道E點附近出現風速增大現象。來流風角為30°和60°時,山體M1和M2對來流風起到了阻擋作用,來流風到達賽道附近觀測點E時風速相應的減小。來流風向角為90°和120°時,來流風經過山體M2和M3出現峽谷風效應,風速在2個山峰之間加速,之后平緩地過渡到觀測點E。來流風向為150°和180°時,來流風經過2個山體,到達E點附近經歷了2次折減。來流風向角為210°、240°、270°和300°時,來流風主要受山峰M4的影響,風向角為270°時阻擋影響較大,觀測點E附近的風速最小。風向角為240°、270°和300°時,來流風均在山體M2和M3之間出現峽谷加速效應。來流風向為330°時,山體M1和M4對來流風起到了阻擋作用,來流風到達賽道附近觀測點E時風速逐漸減小。風的流動方向會受到鄰近山體的影響,在來流風270°較為明顯。整體上看,來流風向為30°、60°、270°和330°時,滑雪場鄰近山體對來流風的阻擋效應最大,賽道處的風速較小,適宜正常開展比賽及訓練。

2.2 索道位置處風環(huán)境評價

索道各位置的風速會受到鄰近山體的影響,從圖15中可以看出,除來流風向為270° 外,來流風速經過復雜地形后,索道各觀測位置的平均風速放大系數Cu基本大于1,風速整體增大?？傮w上,索道各觀測位置風速均受地形地勢條件的影響,地勢較高的觀測位置處的Cu整體大于地勢較低位置處的Cu,但6#觀測點的Cu最小,1#位置的Cu最大。因此,可將1#觀測點位置的風速作為索道安全運營的閾值風速,當此觀測點風速超過索道的安全限值風速20 m/s時[27],應停止索道的運行。來流風向為0°時,索道各觀測點的Cu最大。來流風向為270° 時,各觀測點的Cu最小,分析結果與前文2.1節(jié)中地形對滑雪場區(qū)域風場的影響基本一致。

風攻角是影響纜車氣動力系數的重要參數,也是結構抗風設計中考慮的重要因素[28-29]。由圖16可見風攻角呈180°對稱分布,來流風向為0°～150°時,氣流從低海拔流向高海拔地區(qū),會出現“爬坡”現象,索道1#～9#位置產生正向的風攻角。來流風向為180°～330°時,來流風越過山頂后尾流產生向下的風速分量,氣流出現“下坡”現象,索道1#～9#位置產生負向風攻角。總體上看,各工況索道上站到5#觀測點位置的風攻角較大,這是因為上站到5#觀測點位置之間的海拔較高,地勢與周邊平坦地區(qū)相比起伏大,地形地勢對氣流的影響顯著,產生較大的風攻角。而6#觀測點到索道下站位置之間的海拔相對較低,地勢相對平坦,風攻角較小。在該滑雪場抗風設計中,要考慮風攻角對索道上站到5#觀測點位置之間產生的不利影響。

圖17列出了不同來流風向下索道高度觀測位置的風偏角分布。來流風向為270°時,5#觀測點到索道下站位置之間的風偏角較大,除此來流風向外,風偏角變化較大的位置均在索道上站到3#位置之間,5#～7#位置由于受局部地形的影響,風偏角變化也相對較高。由圖1可知,4#位置由于受到周圍地形地貌的影響各來流風向下風攻角變化均較小,總體上看,索道位置地勢越高對氣流的流動方向影響越大。風偏角對索道纜車的行進安全有重要的影響,根據GB 50127—2020《架空索道工程技術標準》[30]風速大于5 m/s橫向擺動一般小于0.24 rad(約13.75°)時不離開承載索。來流風向為330°時的索道上站到3#位置,以及來流風向為270°時,5#位置到索道下站位置,風偏角均大于13.75°,當風速大于5 m/s時,應考慮風偏角對索道安全運行的影響。

圖17 不同來流風向下索道高度風偏角分布Fig. 17 Wind yaw angle distribution of ropeway height at different incoming wind directions

3 結論

以亞布力滑雪場為工程背景,利用三維地形建模及CFD模擬方法,研究了地形對滑雪場賽道風場的影響以及滑雪場索道區(qū)域風空間分布特性,主要結論如下:

1)文中采用地形過渡曲線對山區(qū)地形進行建模,利用組合網格劃分方式,通過壁面函數方法對研究區(qū)進行CFD數值模擬。數值模擬結果與現場實測的偏差基本在20%以內,CFD數值模擬方法具有一定的準確性,能夠滿足復雜地形滑雪場區(qū)域風特性研究的要求。

2)基于該方法可以有效地開展亞布力滑雪場賽道風場分析,來流風向為30°、60°、270°和330°時,滑雪場鄰近山體對來流風的阻擋效應最大,賽道處的風速較小,適宜正常開展比賽及訓練。

3)該數值方法可以模擬評價亞布力雪場索道風環(huán)境,索道各觀測位置的平均風速放大系數Cu基本大于1(來流風向270°除外)。1#觀測點位置的風速可作為索道安全運營的閾值風速。風攻角和風偏角均受地形地勢的影響,風攻角較大的位置出現在索道上站到5#觀測點之間,風偏角變化較大的位置出現在索道上站到3#觀測點之間。