基于WRF-LES模式的大氣邊界層近地風場精細化模擬研究

劉達琳, 陶 韜, 曹 勇, 周 岱, 韓兆龍

(1. 上海交通大學 船舶海洋與建筑工程學院,上海 200240;2. 安徽工程大學 建筑工程學院,安徽 蕪湖 241060)

近年來,全球性氣候變化使得我國臺風災害頻發,對沿海的房屋建筑造成巨大影響,如臺風“莫拉克”造成的房屋倒塌有1.4萬間[1].臺風等極端氣象對人類產生直接危害的風場來自更加靠近地面的大氣邊界層.研究大氣邊界層,預測邊界層特性,對建筑結構抗風系數設計、土木工程防災減災具有重要意義[2-4].

精細研究大氣邊界層特性常用的風工程方法主要有現場實測、風洞實驗和數值模擬等3類.與現場實測和風洞實驗相比,數值模擬方法不受探測技術、設備性能、地理因素等客觀因素限制,成為使用更多的研究手段.Fluent等經典軟件未考慮地表能量平衡計算、云微物理參數化等氣象模式,因此在真實大氣應用中存在一定缺陷;此外,由于運算能力和時間的限制,Fluent等經典軟件不適合較大范圍的風場模擬,而數值天氣預報(Weather Research and Forecasting,WRF)系統可以彌補這些缺陷.WRF模式具有參數多、方案新、精度高、易嵌套等特點和優點,可用于預測、重現大氣邊界層風場特性.文獻[5]中利用WRF模式數值模擬蘭州河谷盆地區域,獲得的風場和溫度場模擬結果與觀測基本一致,為緩解城市熱島效應提供了參考依據.文獻[6]中采用WRF模式對臺風“鸚鵡”進行高時空分辨率模擬,通過分析臺風模擬路徑、氣壓分布、風速模式等風場特性,發現WRF模式能有效模擬近地面臺風風場,為后續小尺度數值模擬提供真實的入流風場.Yuan等[7]采用WRF模式對真實陸上風電場內的流動特性和功率輸出特性展開高分辨的數值模擬,利用風電場的觀測數據進行驗證,結果表明,WRF模式對于風場的風速、風向和風力機功率的模擬具有較高的精度.WRF模式匯聚了不同學者描述不同物理過程的參數化方案,它們各有特點和使用條件限制,因而會影響數值模擬的穩定性和精度,如空間差分格式方案和次網格模型方案對精細化模擬產生影響.Wicker等[8]對比通量的三階、四階、五階和六階形式,發現奇數階算法是更高一階(偶數階)算法的線性組合,具有數值耗散性.Janjic等[9]在非線性試驗中發現四階差分格式比原始二階差分格式對于小尺度能量積累更具有優勢.Yamaguchi等[10]研究發現高階對流方案有利于減少云模式的數值耗散.文獻[11]中設計25個試驗(水平對流項和垂直對流項分別采用二階～六階5種差分格式),開展非線性密度流試驗研究,通過分析渦旋產生位置和發展,發現水平五階差分和垂直三階差分的組合方案,基本可消除不穩定擾動.大渦模擬(Large-Eddy Simulation, LES)方法通過直接求解大尺度渦,采用次網格模型參數化小尺度渦,有利于精細化模擬近地面風場,Deardoff[12]在1972年將大渦模擬應用于大氣邊界層的模擬.現有WRF版本集合了LES模塊,通過修正WRF-LES模塊的次網格模型而較好地呈現近地面風場特性.Moeng[13]在1984年LES中也采用了TKE閉合模式,但該方法所模擬的平均風速在近地面附近不滿足相似理論的對數律.Moeng等[14]也在傳統的Smagorinsky次網格模型添加修正項,以減小不同計算域間的表面摩擦偏差,改進后的模型較好地重現了熱力驅動或風切變引起的邊界層湍流特性,提高了大渦模擬的精度.Leith[15]考慮了湍能的逆向串級,提出隨機后向散射次網格閉合模式.Mirocha等[16]在WRF-LES中添加非線性回波散射和各向異性(Nonlinear Backscatter and Anisotropy,NBA)模型,與線性模型相比,NBA模型明顯降低近地面風速剖面與對數律預期之間的偏差,從瞬時風場圖中發現其模擬的渦旋尺寸范圍更廣,能更好描述高分辨率下地面附近的流動分離.Kirkil等[17]利用WRF模式研究了拉格朗日平均尺度相關模型和動態重構模型,研究結果表明這兩種動態次網格應力模型能夠很好預測功率譜中慣性區域的產生和范圍擴展,實現最好的動態模型的整體縮放,這種效果在地表附近效果更加明顯.此外,網格分辨率也影響模擬的精度,Zheng等[18]發現提高網格分辨率對空氣污染、降雨降水的預測非常重要.文獻[19]中利用WRF-LES模式評估3種水平網格分辨率(120、60、30 m)和3種垂直網格分辨率(20、10、5 m)對對流層的影響,發現水平網格分辨率為30 m時,能在大氣邊界層范圍獲得理想的湍流再現,而垂直網格分辨率的提高對模擬結果影響較小.Mirocha等[16]發現網格縱橫比(水平網格分辨率/垂直分辨率)為3～4時,平均風速剖面更加接近指數律.

上述研究運用WRF-LES模式模擬大氣邊界層,多采用均勻網格,尚缺少加密網格和非線性次網格模型、空間差分格式奇偶性對近地面風場精細化模擬影響的深入研究,可能會導致該模型難以準確模擬大氣邊界層近地風場特性.

本文基于WRF-LES模式,尋找適用于近地面風場精細化模擬的次網格模型方案、空間差分方案和網格設置方法.首先闡明數值模型的控制方程和湍流模型,然后闡述數值模擬方法,包括大氣邊界層的幾何模型、邊界條件和計算設置,隨后討論和分析了計算結果,最后進行總結.

1 數值模型

1.1 控制方程

利用WRF-LES模式開展理想大氣邊界層近地面風場特性的數值模擬,為了定量描述和預報邊界層狀況,文獻[20]中引入可壓縮流體的狀態方程、連續性方程、動量方程、熱力學方程和水分方程.考慮易讀性,簡化方程表示如下.

狀態方程:

(1)

連續性方程:

(2)

動量方程:

(3)

熱力學方程:

(4)

水分方程:

(5)

式中:ρ為空氣密度;p為壓強;T為溫度;θ為位溫;Rd為比氣常數;xi、xj表示空間坐標,i,j=1, 2, 3;cp為比定壓熱容;ui為流體速度,i=1, 2, 3, 分別表示緯向分量、經向分量以及垂直分量;Sθ為熱量的源匯項,分別由輻射、潛熱、湍流及對流輸送造成;Ωj為地轉角速度在各方向上的分量;Fi為摩擦力項;qn為比濕;Sqn為水分的源匯項,n=1, 2, 3, 分別表示固態水、液態水和氣態水;當i、j、k為相同數值或為偶排列時,εijk為0,當i、j、k為奇排列時,εijk為-1.

將動量方程在空間上進行過濾,只濾去小波脈動,保留大渦脈動,可得:

i,j,k=1, 2, 3

(6)

1.2湍流模型

為了封閉動量方程,WRF4.3提供4種次網格模型:標準Smagorinsky模型(SMAG模型)、標準1.5階湍動能閉合模型(TKE模型)、NBA1模型和NBA2模型.前兩種模型假設次網格應力與渦黏系數是線性關系,而渦黏系數分別由Smagorinsky常數cs和ce來確定;后兩種模型分別是在前兩種模型的基礎上額外增加非線性項,以考慮湍流反向級串和各向異性的影響.NBA1模型和NBA2模型的亞格子應力項[16,21]如下式所示:

(7)

(8)

2 數值模擬方法

2.1 幾何模型與計算網格

通過WRF-LES模式構建模型,以理想大氣邊界層為研究對象,研究次網格模型方案、網格分辨率方案和空間差分格式方案對大氣邊界層近地面風場特性的影響.表1給出理想大氣邊界層的模擬范圍,水平和垂直的網格分辨率和網格數量,以及時間步長.對于次網格模型試驗,水平網格分辨率取30 m,垂直網格分辨率取20 m,選用WRF4.3提供的4種次網格模型.對于網格分辨率試驗,在一個基本算例(水平網格分辨率30 m,垂直網格分辨率 10 m)的基礎上,選用只改變水平網格分辨率的4種試驗方案和只改變垂直網格分辨率的5種試驗方案.其中,不均勻加密方案的計算域如圖1所示,入口風速為15 m/s,計算域長寬為 6 000 m,計算高度為 2 000 m.水平方向上網格分辨率取30 m,網格數為200×200,網格節點均勻分布.垂直方向上采用拉伸網格方案,如圖2所示.圖中:x1、x3分別為經度坐標和垂直坐標.在離地200 m范圍內加密16層,第1層網格高度近似為3.54 m,網格膨脹率為1.1;200 m以上采用間距30 m的均勻網格.垂直網格不均勻加密方法可以在使用少量網格情況下,解決壁面黏性切應力的影響.對于空間差分格式試驗,H代表水平對流,V代表垂直對流,例如H5V3表示為水平對流項取五階差分,垂直對流項取三階差分.

圖1 理想大氣邊界層計算域 Fig.1 Computational domain for ideal atmosphere boundary layer

圖2 垂直網格加密示意圖Fig.2 Illustration of vertical grid mesh encryption

表1 基于WRF-LES的次網格模型方案、網格分辨率分案和空間差分格式的試驗模擬參數Tab.1 Experimental simulation parameters of subfilter-scale stress models, mesh schemes, and spatial difference schemes based on WRF-LES

2.2 計算設置

計算域左側為速度入口,右側為速度出口,前后左右均為周期性邊界條件,頂部邊界條件光滑無滲透,表面層方案采用基于Monin-Obukhov相似理論[22]的改進MM5方案,該相似理論與風工程中常用的壁面函數一致,由以下公式確定:

(9)

(10)

式中:τsurf為地表摩擦力;u*為摩擦速度;ua為z高度處的風速;Cd為摩擦因數;z0為地面粗糙度;κ為von Karman常數,取值0.4;L為Obukhov長度;U為順風向風速;ψ為大氣穩定度的方程[23],在此不贅述.

地面粗糙長度設為0.1 m.時間積分方案采用三階Runge-Kutta方法.所有其他物理選項均關閉,如微物理過程、積云方案、陸面過程、輻射過程等.計算時間為20 h(次網格模型中的NBA1和NBA2方案除外),每隔1 h存儲一次模擬的全流域結果.設立55個監測點,記錄每個時間步監測點剖面的風速,最后2 h的監測點數據用于結果統計分析.將最后2 h的統計數據與更長時間計算得到的數據進行對比,發現二者數據較為吻合,驗證了時長的獨立性.

3 WRF-LES計算方法驗證

根據文獻[16]中的計算設置驗證WRF-LES計算方法的準確性.該試驗中,大氣邊界層高度Ha=1 000 m,水平網格分辨率設為32 m,水平格點數為128,垂直網格分辨率設為8 m,垂直格點數為125,次網格模型采用Smagorinsky方案,時間積分方案采用三階Runge-Kutta,水平對流項采用五階差分,垂直對流項采用三階差分.將計算穩定后的平均風速剖面與文獻[16]中進行對比,如圖3所示.

圖3 試驗模擬和參考文獻[16]結果的平均風速剖面對Fig.3 Validation of the numerical model by comparing the wind speed profiles between test result and reference result[16]

由圖3可知,本節試驗與文獻[16]中的平均風速剖面基本吻合,二者最大的相對誤差為8%,出現在47 m高度處,其余誤差基本小于5%.盡管試驗模擬結果略大于已有文獻,但這種誤差是可以接受的.

4 計算結果和分析

4.1 次網格模型對大氣邊界層模擬結果的影響

圖4給出不同次網格模型下的平均風速剖面.從圖4(a)可以看出,水平風速在近地面呈現近似指數函數的上升規律,大約在120 m高度處逐漸變陡,直到 1 000 m 大氣邊界層高度處停止并再次出現轉折.本文不討論大氣邊界層宏觀空間的問題,而是聚焦近地面的風場特征.近地面層是最接近下墊面的大氣部分,高度通常在50～100 m[24],因此把近地面高度定為100 m.

圖4 不同的次網格模型方案下平均風速剖面Fig.4 Mean wind speed profiles of different subfilter-scale stress models

將圖4(a)平均風速剖面在近地面處進行無量綱化處理,與對數律(log law)進行對比,如圖4(b)所示.在對數坐標下,平均風速剖面與理論曲線吻合較好,初步驗證WRF-LES方法可行性.此外,可看出模擬的水平風速在近地面處均大于理論值;較標準TKE模型或SMAG模型,NBA模擬結果更接近對數律理論值,模擬效果更好.

表2進一步定量分析次網格模型的模擬精度,結果表明NBA1模擬的平均風速與對數律的相對誤差基本最小,與文獻[16]中發現的結果一致.這是因為采用標準TKE模型或SMAG模型時,只能參數化正向的湍流能量串級而無法估算反向的能量輸送,所以會高估近地層的風速切變[25-26].但在NBA模型中考慮次網格隨機擾動作用[27],即考慮了能量可以從小尺度渦旋往大尺度渦旋傳輸;同時引入非均勻因子[28],使得風切變在近地層引起湍流各向異性,導致次網格的湍動能重新分配,從而一定程度上修正風切變對風剖面影響,提高近地風場的模擬精度.

表2 不同的次網格模型方案下平均風速模擬結果與對數律的相對誤差Tab.2 Relative error between the simulated mean wind speed and log law of different subfilter-scale stress models

圖5給出近地面范圍內不同次網格模型的水平方向湍流強度I1,可看出湍流強度基本在0.12～0.22波動,而且30～100 m高度內NBA模型的湍流強度明顯大于標準TKE模型或SMAG模型,在相同的高度處,不同次網格模型間的相對誤差可達到10%,這可能是因為NBA中的非線性成分對來流擾動作用增強,使得湍流強度更大.

圖5 不同的次網格模型方案下湍流強度剖面Fig.5 Turbulence intensity profiles of different subfilter-scale stress models

取所有監測點的順風向風速序列,分別做功率譜,空間平均后再擬合處理,具體結果如圖6所示.圖中:f表示頻率;S(f)為風速譜;σu為順風向速度脈動均方根.從圖中可以看出,頻率處于0.025～0.03 Hz之間,功率譜在雙對數坐標圖中呈現出明顯的線性特性,所有次網格模型都有相似的譜能下降速率,截止頻率都近似在0.03 Hz,但頻率大于0.03 Hz部分,TKE和SMAG線性次網格模型較非線性NBA次網格模型功率譜幅值下降更快.為了進一步定量分析數值解析湍流成分的最小尺度,下文重點關注有效網格分辨率.有效網格分辨率、截斷波長附近的能量衰減模式動能譜顯著偏離實際大氣動能譜的波長.由文獻[27]中分析大氣行為的基本特征,發現在大尺度區域(500～5 000 km)動能譜(E)與波數(r)近似滿足E∝r-3關系,過渡到中尺度區域(400 km以下)近似滿足E∝r-5/3關系,結合泰勒凍結假設,建立功率譜與有效網格分辨率的計算模型如下:

圖6 在100 m高度處,不同的次網格模型方案的功率譜Fig.6 Power spectrum of different subfilter-scale stress models at 100 m above the surface

式中:λe為有效網格分辨率;fe為有效頻率.

根據上式計算得有效網格分辨率近似為12Δx.但在其他文獻中,文獻[29]中采用動能譜評估WRF模式,確定7Δx為WRF模式的有效網格分辨率;文獻[30]中也發現WRF-LES模式可以分辨7倍格距波長.但本文得到的有效網格分辨率較高,可能是因為本文通過功率譜與斜率為 -5/3 的直線嚴格相切來確定有效網格分辨率,而文獻[27-28]中則依據斜率變化大致進行判斷.另外參考文獻[30-31],還對水平分辨率為8.2、30、90 和270 m的功率譜進行嚴格相切,發現有效網格分辨率近似在11～15倍格距,與本文結果基本一致.

因此,綜合水平風速剖面、湍流強度剖面和有效網格分辨率等方面考慮,確定NBA1模型為次網格模型的最佳選擇.

4.2 網格分辨率對大氣邊界層模擬結果的影響

網格分辨率對數值模擬的影響是一個相對復雜的問題.一方面數值模擬精度依賴于網格分辨率,為提高數值模擬精度要求網格分辨率盡可能高,尤其是在壁面附近;另一方面,受計算條件限制,網格分辨率不能無限增加.因此,討論如何在節約計算資源下優化網格設計,使得模擬結果真實可靠.

圖7 不同水平網格分辨率和垂直網格分辨方案的平均風速剖面Fig.7 Mean wind speed profiles of different mesh resolution cases

表3 不同水平網格分辨率和垂直網格分辨率方案下平均風速模擬結果與對數律的相對誤差Tab.3 Relative error between the simulated mean wind speed and log law of different horizontal and vertical mesh cases

由圖7(b)可看出,盡管垂直網格不均勻加密對平均風速剖面影響不大,但在近地面處會使其結果更靠近對數律,結合表3,當垂直網格最小尺寸越小,近地面平均風速的相對誤差會一定程度上減小.此外,當垂直網格分辨率從 30 m 加密到 20 m (網格數為400萬)時,模擬結果已經開始收斂,若采用垂直網格不均勻設計(網格數為308萬),網格總量降低了23%,顯著降低計算規模,在滿足精度要求下提高計算效率,因此有必要加密近地面網格.并且當網格縱橫比為3,相對誤差能控制在10%以內,模擬效果較好.

圖8給出不同網格分辨率下的湍流強度剖面,可以看到湍流強度基本在0.12～0.24波動,但在不同高度處湍流強度變化存在差異;當采用相同的垂直網格分辨率時,如圖8(a)所示,湍流強度在一定高度范圍內隨水平網格分辨率增加而增加,湍流強度拐點高度近似都為15 m;當采用相同的水平網格分辨率時,如圖8(b)所示,湍流強度在一定高度范圍內也隨垂直網格分辨率增加而增加.而超出某個高度,湍流強度會隨水平網格分辨率或垂直網格分辨率增加而減小,這可能因為當網格分辨率越高,模擬的平均風速越接近入口風速,即平均風速越大,從而湍流強度減小.此外,如圖8(b)所示,湍流強度拐點高度會隨垂直網格分辨率變化而變化,當垂直網格分辨率為10 m時,湍流強度拐點高度近似為5 m,當垂直網格分辨率為30 m時,湍流強度拐點高度近似為 45 m;模擬結果的湍流度剖面拐點通常出現在地面以上第2層網格處.原因在于真實大氣的湍流度越靠近地面越大,然而在模擬中,由于靠近地面網格分辨率和壁面邊界條件的限制,第1層的湍流難以充分直接解析,所以最高湍流度往往出現在第2層.

圖8 不同網格分辨率方案下的湍流強度剖面Fig.8 Turbulence intensity profiles of different mesh resolution cases

綜上考慮,推薦水平網格分辨率為30 m,垂直網格分辨率為10 m,垂直網格在近地面區域適度加密.

4.3 空間差分格式對大氣邊界層模擬結果的影響

空間差分格式會影響N-S方程對流項的求解,為了解空間差分格式對理想大氣邊界層數值模擬影響,圖9給出不同空間差分格式的平均風速剖面.可知,所有空間差分格式的無量綱化平均風速基本都在14.5～18.0波動,與對數律非常接近,所有相對誤差基本控制在5%以內,說明空間差分格式對平均風速剖面影響不大.

圖9 不同空間差分格式方案下的平均風速剖面Fig.9 Mean wind speed profiles of different advection differential cases

圖10給出不同空間差分格式在近地面范圍的湍流強度剖面.由圖可知,空間差分格式在約40 m 高度范圍內存在明顯的差異,奇數階差分格式較偶數階的湍流強度更大;但在 40 m 高度以上,所有空間差分格式的湍流強度都趨于0.16.

圖10 不同空間差分格式方案下的湍流強度剖面Fig.10 Turbulence intensity profiles of different advection differential cases

圖11給出不同空間差分格式在同一時刻的水平速度,剖面高度為z=40 m.由圖可見,所有試驗方案的流場圖均有條帶狀結構,湍流結構有明顯的組織性,而且受地轉風影響,條帶狀結構呈現近似45° 方向.另外,由圖可看到,偶數階差分格式較奇數階明顯能捕獲更小尺度的湍流結構,出現這種差異的原因可能是奇數階差分格式較偶數階增加了數值黏性,使得小尺度渦耗散,小尺度湍流結構消失.結合圖12和表4,可得偶數階差分格式的有效分辨率為(6～8)Δx,奇數階的有效分辨率為(9～13)Δx,證實偶數階差分格式能獲得更小尺度的湍流結構.

圖11 在40 m高度處,不同空間差分格式方案下x-y平面水平速度的瞬時風場Fig.11 Contours of instantaneous velocity in the x-y plane at 40 m above the surface from different advection differential cases

圖12 在40 m高度處,不同空間差分格式方案下的功率譜Fig.12 Power spectrum of different advection differential cases at 40 m

表4 不同空間差分格式方案下,40 m高度處的有效網格分辨率Tab.4 Effective resolutions of the different advection differential cases at 40 m

5 結論

利用WRF-LES模式對理想大氣邊界層進行模擬,結合平均風速剖面、湍流強度和功率譜等數據,分析次網格模型、網格分辨率和空間差分格式對近地面風場模擬的影響.主要結論如下:

(1) 與標準TKE模型或SMAG模型相比,NBA1次網格模型考慮了湍動能的反向級串,可更充分揭示湍流性態,更準確地模擬理想大氣邊界層近地面的風場特性.

(2) 結合水平網格分辨率(15、30、60和120 m)和垂直網格分辨率(5、10、20、30 m和不均勻加密)試驗方案,評估其在近地面的模擬精度,建議網格縱橫比設為3,水平網格分辨率采用30 m,垂直網格分辨率采用10 m;垂直網格在近地面區域適當加密可有效節約計算資源.

(3) 空間差分格式對平均風速剖面影響不大,但在WRF-LES模式下,與奇數階差分格式相比,偶數階差分格式,可明顯改善小尺度湍流結構的再現能力,即偶數階差分格式有效分辨率達(6～8)Δx,而奇數階差分格式的有效分辨率為(9～13)Δx.

合適的WRF-LES參數方案可顯著提高近地面風場模擬精度,有望應用于臺風邊界層的精細化數值模擬,提高臺風風速預測精度,為防臺風減災設計提供技術參考.