一種山區峽谷橋址區風場特性數值模擬方法

靖洪淼， 廖海黎,2， 周 強,2， 馬存明,2

(1.西南交通大學 土木工程學院，成都 610031；2.風工程四川省重點實驗室，成都 610031)

近年來，山區建造大跨度橋梁、大跨空間結構等風敏感性基礎設施的情況逐漸增多。對于此類柔性結構而言，建設場地的強風特性是決定其風致響應的關鍵因素之一，有時甚至是控制其結構設計的決定性因素。現有的相關抗風設計規范[1]也大多針對平原、緩丘和沿海地區，不能直接用于確定復雜山區峽谷地形橋址處的風場特性。對此，目前主要有三種方法用于研究并確定復雜地形的風場特性，分別是：現場實測、地形模型風洞試驗和CFD(Computational Fluid Dynamics)地形數值模擬。

現場實測方法是直接最有效獲取風場特性參數的手段，尤其是脈動風速特性中的風速譜、紊流強度、紊流積分尺度、陣風系數等。例如張明金[2-4]等，在現場建立測風塔，安裝風場實時監測系統，得到了較為完備的山區峽谷地形橋位處的風場特性。但是，該方法需要很長的數據采集時間，現場條件苛刻，耗費大量的人力、物力和財力。地形模型風洞試驗很好的避免了現場實測方法的缺點，如試驗周期短，試驗條件容易控制，模型制作簡便，試驗數據可靠等。Xu等[5-9]采用地形模型風洞試驗方法，研究了復雜山區峽谷地形橋址處的風場空間分布特性。受風洞試驗條件的限制，地形模型表達的地理范圍和采用的比例尺一般都很小。此外，地形模型風洞試驗無法保證雷諾數一致性要求，其試驗結果準確性還有待進一步研究。CFD地形數值模擬很好地解決了現場實測和地形風洞試驗中的缺點，如地形范圍可任意選取，風速和風向可任意設置，費用最少，周期短等。目前采用CFD方法研究風場特性也是一大熱點，例如李永樂等[12-13]采用了地形CFD數值模擬方法，得到了較大范圍復雜地形的風場特性。Ren等[14]利用數值模擬方法，研究了河南某地的風場特性，并在此基礎上提出了一種新的復雜地形空間風場預測模型。Jubayer等[15]同樣采用數值模擬方法，研究了一個較大區域地形風場特性，并以數值模擬結果為基礎設置地形風洞試驗入口條件，最終得到了更加接近實際情況的風場特性。此外，部分文章還對地形過渡段形式做了研究并應用，例如胡朋等[16]基于理想流體圓柱繞流推導出一類過渡段曲線， Ren等采用了一種平方余弦函數用于計算域邊界平滑過渡。

雖然CFD地形數值模擬得到了普遍的研究和應用，但入口風速條件設置、地形計算區域選取、地形邊界過渡形式等方面還不明確，如入口的參考風速仍然很隨意，使得數值模擬雷諾數和實際雷諾數嚴重不相符；截取矩形地形作為計算區域，難于保證風場的充分發展，這些都會導致計算結果的合理性和準確性不足。為此，本文提出了一種新型的山區峽谷橋址處風場特性數值模擬研究方法，合理地設置了計算域，確保風場充分發展；并采用“虛擬標準氣象站”法獲取較合理的入口邊界條件，使數值模擬雷諾數與實際雷諾數盡量一致。最終，根據本文提出的數值模擬方法，得到了較為準確的山區峽谷橋位處風場特性參數，也為該類地形數值模擬研究提供了一定參考。

1 數值計算模型

某待建大跨度橋梁走向以及橋位附近的地形如圖1所示，從圖中可以看出，河道總體呈西南到東北走向，在橋位處幾乎與待建橋梁垂直。橋位四周均有海拔達1 000 m以上的山峰，并且數量眾多。橋位附近東部和北部相對平坦，海拔較低，而遠處卻是海拔較高山峰。因此該地區地形十分復雜，其風環境也非常復雜。

1.1 地形和計算域

本文采用imageware軟件，直接利用散點云圖數據生成地形曲面，即首先基于SRTM(Shuttle Radar Topography Mission)數據(SRTM數據是根據美國太空總署、國防部國家測繪局以及德國與意大利航天機構共同合作完成聯合測量的雷達影像數據，繪制成的數字地形高程模型)，在專業地形軟件Global Mapper中，提取三維地形散點云圖，然后，導入imageware軟件，直接利用上述散點云圖構建三角形網格地形曲面。最后在地形曲面上選取以橋位為中心的半徑5 km的圓形地形范圍。其中三維地形散點云圖的水平解析度是30 m，以及圓形地形范圍內最高和最低海拔相差約750 m，imageware軟件擬合的地形曲面與原數字高程數據的誤差不超過5%。

圖1 橋位處地形圖Fig.1 Terrain at bridge site

由于橋位處地形非常復雜，而且地形解析度有限，所以導致計算域邊界十分陡峭。為減少地形邊界突變對風場模擬結果的影響，確保獲得更符合實際的風場選取的地形邊界需要做延伸以及光滑處理。本文利用ICEM軟件中的B-Spline曲面形式作為過渡面，連接地形邊界與計算域底面，如圖2所示。此外由于選取的地形區域面積足夠大，當氣流發展到橋位處時，邊界過渡段的影響基本可以忽略。在本研究中，地形計算域的過渡段的距離為1 km，而Liu 等[17]在30 km×30 km的計算域中僅采用了1 km長度的過渡段，并取得了較好的模擬結果。顯然本文計算域遠遠小于Liu等采用的計算域，但仍然采用1 km的過渡段，因此上述處理方式可以滿足數值模擬的需要。最終處理后的地形如圖3所示，其中待建橋梁主梁和地形邊緣的海拔高度差為70 m。

為了確保風場的充分發展，并兼顧工況設置，總體計算域在圓形地形區域上，采用長方體形式，長寬高分別為：30 km，20 km和3 km，如圖4所示。計算域的阻塞率不超過5%。所關注的地形區域(半徑為6 km的圓形區域)位于總長方體計算域中心，而且可以繞圓心旋轉，便于改變來流方向，從而滿足設置不同計算工況的要求。

圖2 地形邊界過渡形式Fig.2 Transition of terrain edge

圖3 地形三維曲圖Fig.3 Three-dimensional terrain

圖4 計算域Fig.4 Computational domain

1.2 網格劃分

如圖5所示，整個計算域分成兩個部分：外部區域和內部圓柱地形區域。其中外部區域比較規則，因此采用正交性較好的六面體結構化網格劃分；內部地形很復雜，因此采用適應性更好的四面體非結構化網格劃分。內部地形表面采用了10層邊界層網格，其中第一層網格的水平間距不超過40 m。針對第一層網格豎向高度，根據相關文獻給出的經驗值(如李永樂等、Ren等、Liu等和祝志文等[18]給出的第一層網格高度分別為10 m，5 m，6 m和5 m)，本文設定為5 m，故其y+值為500左右。Blocken[19]指出y+≈500以內屬于邊界層，且在邊界層內可采用壁面函數處理，節省計算資源而且獲得較精確的計算結果。因此，本文在地形表面增設了壁面函數(nut-K wall function)，以確保數值計算的精度。此外，Liu等指出地表樹木等植被的平均高度為15 m，第一層網格高度小于此值即可滿足要求，即本文設置的第一層網格高度滿足此要求。由于地形復雜，在內部圓柱地形區域采用非結構化的四面體網格，因此，地表附近部分網格質量稍差，但質量也在0.35以上，其余部分網格質量均大于0.5。外部區域采用結構化六面體網格劃分，底層網格高度為5 m，從下往上也是按照1.2的比例逐漸增大。總體網格數量約為420萬。

1.3 湍流模型和邊界條件

在風場數值模擬中，選取了對流動分離具有較好解析度的k-ωSST湍流模型。采用有限體積法對計算域流場進行離散，由于復雜地形下橋梁抗風安全性更關心平均風特性，而且計算流場充分發展后的平均風場特性接近穩態計算結果，因此本文采用了穩態不可壓縮simpleFoam求解器。壓力與速度耦合采用simple算法。對流項和黏度項采用了二階中心差分格式離散。

如圖4所示，上面、下面和頂面為對稱邊界條件(symmetry),地面采用無滑移固定壁面(No-Slip)。由于本計算采用的計算域較大，在出口邊界上流場已得到充分發展，故出口邊界采用開放式零壓力出口(Open pressure)。由于橋梁主梁距地表的高度遠遠大于地表植被等的高度，因此在數值模擬中并未考慮地表粗糙度。內部和外部網格交界面采用cyclicAMI邊界條件，該種邊界條件允許速度和壓力的傳遞。計算域入口為速度入口(Velocity inlet)邊界條件，按大氣邊界層C類指數風剖面(粗糙度系數為0.22)形式從計算域入口底面開始設置；在梯度風高度450 m以上，風速采用固定值。

為了更加科學地設置入口橋位高度處的參考風速，使模擬地形風場的雷諾數更加符合實際情況，采用了“虛擬標準氣象站”法[20]計算入口橋位高度處的風速值。現假定在橋位主梁高度處設立虛擬氣象站，該站的平均海拔高度取為688 m。從《公路橋梁抗風設計規范》中，可以查出橋位200 km范圍內8個地區的海拔和對應的百年風速，見表1。圖6給出了表1中各氣象站點最大風速隨高度變化的情況，并利用最小二乘法對數據進行擬合，得到橋位處虛擬氣象站重現期為100年的風速隨海拔高度變化的曲線。利用擬合的風速-海拔高度公式，可得該虛擬氣象站的百年一遇的最大風速為24.75 m/s。將24.75 m/s作為橋位虛擬氣象站的基本風速U10。入口采用指數風速形式，根據橋位周圍的地形情況，取為C類地形，即地面粗糙度系數α=0.22，湍流度取5%。入口指數風速表達式為式(1)，如前面所述，橋梁主梁到計算域入口底邊的相對高度是70 m，所以入口橋位高度處的風速為37.97 m/s。

(1)

式中：z為距計算域入口底邊的距離,m;Uz為距計算域入口底邊距離z處的風速,m/s。

圖5 計算域網格信息Fig.5 Mesh system

地點距離/km海拔/m100年風速貴陽1611 074.325.2桐梓12497225.1習水1881 180.225.4遵義市98843.924.9思南54416.324.4銅仁130279.724.2凱里市109720.324.8吉首市184206.624.2

1.4 工況說明

在實際情況中，影響橋梁工程安全的主要是橫橋向來流，同時為了說明本文數值模擬的合理性，這里設置了6種不同來流風向，見圖7所示，并以此設置工況編號(圖中圓圈內數字)。工況1、工況2和工況3的水平夾角依次為15°；工況4、工況5和工況6的水平夾角也依次為15°。主梁考察點見圖8所示，共18個，依次均勻分布在兩個橋塔之間的主梁上。兩個橋塔間距為680 m，因此主梁相鄰觀測點間距為40 m。此外，為了得到橋位處風速剖面情況，分別在1/4跨、跨中和3/4跨處從下往上依次設置了30個觀測點。其中相鄰觀測點之間的豎向間距為20 m，最底部觀測點到主梁的距離為60 m。

圖6 橋位附近氣象站點最大風速隨高度的變化Fig.6 Variation of maximum wind speed with height of weather stations near the bridge site

圖7 工況示意圖Fig.7 Diagram of load cases

圖8 主梁觀測點布置示意圖Fig.8 Diagram of measurement points at the bridge girder

2 結果與討論

由于本文采用的是穩態計算，故選取流動充分發展的狀態為最終流場狀態。如圖9所示，中心平面上遠離地形的尾流區域的壓力梯度約為零，即表明流動已得到充分發展。

圖9 中心剖面壓力云圖Fig.9 Pressure contour of center plane

基于上述最終流場狀態，分析了山區峽谷橋址處的風場特性，得到了擬建橋梁主梁上的水平風速、風攻角α、風向角β和風速放大系數Cu等，以及主梁1/4跨、中跨和3/4跨處的平均風速剖面情況，并結合CFD數值模擬流跡顯示的優點，研究了橋址處山體等地形對風場特性的影響機理。風攻角α、風向角β和風速放大系數Cu的計算表達式見式(2)～式(4)。

(2)

(3)

(4)

式中：u，v，w分別為主梁上橫橋向，順橋向和豎向風速分量；vb為主梁上的測點風速，v0為入口主梁高度處的風速。

2.1 平均風速沿主梁分布

為了便于橋梁工程抗風設計的需要，主梁水平方向的風速分解為順橋向和橫橋向，并無量綱化，參考風速為主梁高度入口處的風速，如圖10和圖11所示。從圖10的主梁橫橋向風速可以看出，工況4和工況6中部分橫橋向風速較大，其中工況4為主梁前半段，工況6為主梁后半段，另外工況5的橫橋向風速普遍較大。以工況5為例，從圖12的主梁流線圖中可以看出，由于河道和兩岸地勢的原因，來流在主梁處發生了匯聚，產生了明顯的“峽谷效應”，工況4和工況6的主梁流線圖也基本同上。另外，工況2的主梁左半段風速明顯相對較小，分析其主梁流線圖(見圖13)可以發現，由于其上游地形比較復雜，有兩股來流匯聚到主梁上，其中主梁左半段由上部來流控制。這股來流在上游發生了翻轉，特別是到達主梁前，較高的地勢對來流產生了阻擋和分散作用，因此導致風速變小。而下部來流在河道均處于匯聚狀態，流經主梁時處于山體迎風坡頂部，產生了加速作用。

圖10 主梁橫橋向風速分布Fig.10 Distribution of transverse wind speed along the bridge girder

圖11 主梁順橋向風速分布Fig.11 Distribution of longitudinal wind speed along the bridge girder

圖12 工況5的主梁流線圖(左側為風速入口)Fig.12 Streamline of case 5 (inlet on the left)

對于順橋向風速，隨著來流和橋梁橫橋向夾角的增大，風速也由小變大，并且整體分布比較均勻。由于地形相對平緩，該段河道相對平直，因此上下游來流時，順橋向風速線型比較一致，大體呈對稱分布，如工況1和工況6。

圖13 工況2的主梁流線圖(左側為風速入口)Fig.13 Streamline of case 2 (inlet on the left)

2.2 平均風速沿高度分布

為了研究山區峽谷橋址處平均風速沿高度的分布情況，圖14、圖15和圖16分別列出了主梁1/4跨、中跨和3/4跨處，從主梁下60 m到主梁上520 m的風速剖面。從圖中可以看出，當河道下游來流時(對應的工況為4、工況5和工況6)，上述三個跨度的風剖面都比較符合指數形式。這是因為下游地形相對海拔較低，地勢相對平坦，沒有突出的山峰或峭壁，河道也比較平直，因此整個地形對來流的阻擋作用很小，最終導致橋位處的平均風速風剖面形式比較符合指數形式。但是，當河道上游來流時(對應的工況為1、工況2和工況3)，上述三個跨度的風剖面變化非常大，有些已經完全無法用指數風剖面形式表示。產生這種情況的原因是，河道上游地形十分復雜，并出現了河道轉折，在轉折處又存在懸崖，河谷、山谷和峭壁交錯出現，使來流發生翻轉、扭曲、匯聚、擠壓和越山加速等情況，如圖13所示。對于工況1和工況3，由于來流經過較高的山體時發生了山頂阻擋現象，而主梁剛好位于山體背風側，因此部分位置的平均風剖面在低海拔處出現了減速情況。以工況1中主梁3/4跨處的平均風剖面為例，從圖17顯示的主梁流線圖上可以清楚的看到，在主梁的3/4跨附近(圖中箭頭一側)，來流受到河岸突出的山體阻擋而發生跨越，當來流到達主梁時剛好處于山頂背風側附近，因此風速發生了明顯的減速作用。

圖14 主梁1/4跨處風速沿高度變化Fig.14 Variation of wind speed with height at 1/4 span

隨著高度的增加，地形的影響因素也漸漸地變小，平均風速總體也呈現出隨高度增加而增大的現象，并趨向于指數律分布形式，這也說明數值模擬的計算域設置合理，計算域頂部邊界對內部流場的影響很小。

圖15 主梁跨中風速沿高度變化Fig.15 Variation of wind speed with height at mid span

圖16 主梁3/4跨處風速沿高度變化Fig.16 Variation of wind speed with height at 3/4 span

圖17 工況1的主梁流線圖(左側為風速入口)Fig.17 Streamline of case 1 (inlet on the left)

2.3 主梁風攻角和風向角分布

風攻角是影響主梁截面氣動力的重要參數，也是橋梁抗風設計的重要依據。圖18為不同工況下主梁上的風攻角分布情況。總體上看，橋梁主梁的前半段風攻角較小，后半段風攻角較大。這是因為前半段的海拔較低，地勢相對平緩，后半段的海拔較高且地勢起伏大，并且橋梁主梁距離地面相對較近，因此當氣流經過時會受到較為顯著的地形影響，產生較大的風攻角。當河道上游來流時，工況1、工況2和工況3的正負風攻角相對較大，正攻角最大約+10°，負攻角最大約-9°，說明氣流經過主梁時，會分別產生向下和向上的氣流。從圖19中工況3的主梁流線圖上看到，主梁流線有兩股匯聚而成。受上游山體分割，上半部分來流經過河道到達主梁，受河道地勢的干擾，氣流產生向上的流動趨勢，因此產生了正攻角；下半部分來流越過上游山體到達主梁，下半段主梁剛好處在背風坡，因此產生了負攻角。而下游來流時，工況4、工況5和工況6均會產生較大的正攻角，說明氣流在主梁處產生向上流動的趨勢。從地形上看，氣流從地勢低處流向高處，主梁處在迎風坡，自然產生正風攻角。總之，山區峽谷橋址處的主梁風攻角變化非常復雜，局部地形對其有主導作用，因此在實際橋梁抗風設計時，還要結合具體橋址處周圍的地形而確定，特別是當橋梁位于較為陡峭的地形時。

風向角也對橋梁安全有重要的影響，特別是橋梁建設階段。圖20列出了6個不同工況下的主梁風向角分布情況。由于地形高差相對不大，雖然有山體對氣流產生阻擋作用，但距離橋位較遠，氣流繞過山體到達主梁時方向趨向一致，因此導致主梁處氣流方向基本和來流方向一致，除了靠近右側陡峭山體部分外，主梁上的風向角變化不是很大。因此，在該橋梁抗風設計時，需重點關注主梁右半部分的風向角影響，而左半部分按和來流方向一致處理即可。

圖18 主梁風攻角分布Fig.18 Distribution of wind attack angle

圖19 工況3的主梁流線圖(左側為風速入口)Fig.19 Streamline of case 3 (inlet on the left)

圖20 主梁風向角分布Fig.20 Distribution of wind directional angle

2.4 風速放大系數

為了研究該橋址處河道峽谷對氣流的加速作用，按式(4)定義了風速放大系數。圖21為不同來流工況下，主梁上的水平風速放大系數分布情況。從圖中可以發現，工況1和工況6的風速放大系數最大，其中最大值為1.303 2。其它工況的風速放大系數也較大，其中工況2的風速放大系數最小，在2.1節中分析了原因。由于該橋位處的河道較為平直，一定范圍內河道幾乎沒有彎曲，同時地形起伏相對較緩和，海拔差距不是很大，所以不同工況下的風速放大系數較為一致。受河道和兩岸山體的作用，氣流發生匯聚加速，在主梁處產生了較大的風速放大系數。

圖21 主梁風速放大系數分布Fig.21 Distribution of speed-up factor

為了對比說明本文提出的數值模擬方法的正確性，另外采用了規范中的經驗關系式(5)計算橋址處的山谷效應修正系數，其中假設大橋橋址“虛擬標準氣象站”的基本風速為峽谷進口風速。由于該橋址河道較為規則，因此較適用經驗公式計算峽谷效應修正系數。

(5)

式中：H為峽谷深度，一般橋梁建于峽谷上口處，可取主梁橋面至峽谷下底面高度；B1為峽谷上口處寬度，一般為主梁橋面長度；B2為峽谷下底面寬度；κ為山谷效應修正系數。

由于河道分為有水和無水兩種情況，而本文CFD數值模擬計算的是無水情況，因此下面計算河道無水情況。根據式(5)，其B1=680 m，B2=0，H=310 m，π=3.141 6，可得κ=1.286 4。CFD數值模擬得到的風速放大系數最大值為1.303 2，可得與峽谷效應修正系數κ=1.286 4的誤差僅為1.3%，在一定程度上說明了本文提出的數值模擬方法的合理性和準確性。

3 結 論

本文以某在建的山區峽谷大橋橋址處風場為研究對象，提出了一種改進的山區峽谷地形風場數值模擬方法。該方法的計算域設計更為合理；地形邊緣過渡段較為平緩，氣流受到干擾較小，更符合實際流動情況；同時，利用“虛擬標準氣象站”法獲得了計算域入口風速條件，使得計算結果更加準確。通過分析主梁不同位置上的橫橋向和順橋向風速及其放大系數、風攻角、風向角和平均風速剖面等風場關鍵參數及其流動機理，得到了以下結論：

(1)本文采用的改進的山區峽谷地形風場特性CFD數值模擬方法，利用“虛擬標準氣象站”法計算橋址處基本風速，并以此設計計算域入口指數風速剖面，同時設置地形過渡段等改進措施，是正確的、合理的。

(2)采用“虛擬標準氣象站”法計算得到的橋址處基本風速為24.75 m/s。數值模擬結果顯示，主梁上的橫橋向風速沿橋跨變化較大，順橋向風速沿橋跨變化較小；受山地地形的影響，風剖面已經不符合規范中規定的指數風剖面；該橋址處的主梁風攻角在-9°～+10°。

(3)數值模擬得到的主梁上的風速放大系數最大為1.303 2，同時根據規范中的經驗公式得到的山谷效應修正系數為1.286 4，兩者吻合較好，誤差僅為1.3%。