基于WRF與CFD的復(fù)雜地形風(fēng)場 多尺度耦合模擬分析

楊易 譚健成 金博崇 劉慕廣

(華南理工大學(xué) 亞熱帶建筑科學(xué)國家重點實驗室，廣東 廣州 510640)

隨著我國社會經(jīng)濟和建筑行業(yè)的快速發(fā)展，高層建筑的數(shù)量和高度不斷增加。一方面，風(fēng)荷載涉及建筑結(jié)構(gòu)抗風(fēng)的安全性問題，例如近年來不少沿海地區(qū)的高層建筑，在強風(fēng)作用下出現(xiàn)風(fēng)振加速度不能滿足舒適性要求、圍護結(jié)構(gòu)遭到不同程度破壞的問題；另一方面，風(fēng)荷載取值也關(guān)系到建筑結(jié)構(gòu)設(shè)計和建造的經(jīng)濟性問題。因此，作為建筑結(jié)構(gòu)風(fēng)荷載的輸入源頭，大氣邊界層自然風(fēng)場特性的研究對高層建筑結(jié)構(gòu)抗風(fēng)的精細化分析越來越重要。邊界層風(fēng)場是一個涵蓋了幾米到幾千米的多尺度問題，我國現(xiàn)行建筑結(jié)構(gòu)荷載規(guī)范給出了4類不同粗糙度地貌風(fēng)場特性定義，對于平均風(fēng)剖面，梯度風(fēng)高度以下采用指數(shù)率模型描述，但真實復(fù)雜地形的風(fēng)場特性是如何分布，梯度風(fēng)高度以上風(fēng)場又是怎樣的規(guī)律，還亟待研究。

自20世紀80年代以來，中尺度大氣數(shù)值模式已經(jīng)有了很大的發(fā)展，整合了各種最新研究結(jié)果開發(fā)的WRF(中尺度天氣預(yù)報)模式因其開源、高效、可擴展等優(yōu)點已得到了越來越廣泛的應(yīng)用，已有研究表明WRF能夠較為準確地模擬得到臺風(fēng)路徑和復(fù)雜地形風(fēng)場結(jié)構(gòu)特性[1- 4]。黃銘楓等[1]基于WRF與LES(大渦模擬)對九堡大橋周邊臺風(fēng)風(fēng)場進行了多尺度耦合模擬，模擬得到的瞬態(tài)風(fēng)速結(jié)果與實測數(shù)據(jù)接近；劉郁玨等[4]對復(fù)雜山地地形區(qū)域進行WRF-LES在線耦合模擬，并匹配高分辨率的地形高程數(shù)據(jù)，顯著提高了風(fēng)結(jié)構(gòu)的模擬效果，更好地捕捉小尺度風(fēng)速脈動特征。

四維數(shù)據(jù)同化方法(FDDA)[5- 6]，也稱松弛方法(Nudging)，是中尺度數(shù)值模擬的常用加載模塊，它能有效減小WRF模擬過程中的累積誤差。四維數(shù)據(jù)同化方法主要分為以下4類：分析松弛法(AN)、近地面分析松弛法(SAN)、譜松弛法(SN)以及觀測數(shù)據(jù)松弛法(ON)[7]。Tse等[8]成功通過加載FDDA模式提高了WRF模擬結(jié)果的準確性，使得臺風(fēng)邊界層風(fēng)速剖面與實測數(shù)據(jù)有較好的一致性。González-Alonso de Linaje等[9]對比了SN方法和未加載四維同化模塊的控制試驗，結(jié)果表明同化后的均方根誤差和判定系數(shù)都有一定程度的改善。

中尺度氣象模式WRF的模擬區(qū)域大(千米量級)，還考慮了水汽、長短波輻射、積云和下墊面等物理過程的影響，而CFD則具有時空分辨率高、能建立高精度的真實建筑和地形模型的優(yōu)點[10- 11]。國內(nèi)外已有不少研究通過區(qū)域嵌套的方法耦合兩者，綜合兩者的優(yōu)勢，獲得了更準確的精細化風(fēng)場信息[12- 16]。Han等[14]利用WRF輸出地區(qū)平均風(fēng)速的入口邊界值，然后采用BPI(塊多項插值)方法將WRF數(shù)據(jù)插值到LES計算域入口邊界上，模擬結(jié)果表明數(shù)值模擬預(yù)測與實測數(shù)據(jù)吻合良好。Mughal等[15]基于肯尼亞北部圖爾卡納湖所處復(fù)雜地形，通過調(diào)試不同的物理方案和初始場進行風(fēng)場模擬研究，結(jié)果表明，通過WRF-LES模式輸出的風(fēng)速時程結(jié)果相較WRF模擬提高了10%。

本文選取高層建筑密集的深圳地區(qū)作為研究對象，基于中尺度模式WRF與CFD的耦合模式對深圳氣象塔鄰近復(fù)雜地形區(qū)域進行風(fēng)場模擬研究(由于可以和氣象塔風(fēng)場實測數(shù)據(jù)進行對比)，對比分析加載不同同化方案的WRF模擬結(jié)果，依靠中尺度WRF模擬所得環(huán)流信息為小尺度CFD風(fēng)場提供真實的邊界條件；同時加載了更高精度的地形高程數(shù)據(jù)，從而實現(xiàn)降尺度精細化風(fēng)場模擬；最后通過與某一典型風(fēng)場事件下深圳氣象塔觀測數(shù)據(jù)的對比，對數(shù)值模擬的準確性與適用性進行評價。

1 WRF中尺度風(fēng)場模擬

1.1 WRF模式設(shè)置

本文選取北京時間2019年9月12日20時至2019年9月14日2時這一時間段，作為WRF模擬的起止時間(由于此時間段內(nèi)研究場地的風(fēng)場條件相對較穩(wěn)定)，以深圳氣象觀測塔位置(113.898 3°E，22.646 8°N)作為模擬的區(qū)域中心，WRF模擬區(qū)域及氣象塔附近的地貌特征如圖1所示。氣象塔共設(shè)置了13個觀測平臺(具體高度H見表1)，可獲得風(fēng)速、風(fēng)向等實測數(shù)據(jù)。

(a)WRF模擬區(qū)域劃分圖

(b)d03區(qū)域地形高程圖

(c)谷歌衛(wèi)星地圖

(d)氣象塔全貌

(e)氣象塔周邊地貌

(f)氣象塔平臺及伸臂架設(shè)

表1 氣象塔平臺布設(shè)高度Table 1 Height of the platforms at the tower

WRF模式采用3層雙向嵌套的網(wǎng)格劃分方案(d01、d02、d03)，嵌套網(wǎng)格的具體參數(shù)設(shè)置見表2，dx為x方向網(wǎng)格距離，dy為y方向網(wǎng)格距離；網(wǎng)格垂直方向則劃分為50層，貼近地面的20層采用加密網(wǎng)格，首層網(wǎng)格高度約為離地10 m。選用每6 h 更新一次的NCEP 0.25°×0.25°的氣象資料[16- 17]作為模式的初始場和邊界場，地形資料由30 s分辨率的GMTED2010數(shù)據(jù)(The Global Multi resolution Terrain Elevation Data 2010)提供，圖1(b)為d03區(qū)域地形高程示意圖；地圖投影采用Mercator方案。d01、d02區(qū)域結(jié)果輸出的時間間隔均取 1 h，d03區(qū)域則取6 s。根據(jù)過往相關(guān)學(xué)者的研究經(jīng)驗和建議，WRF模擬的具體物理方案設(shè)置如下：微物理過程選用WSM3方案，長波及短波輻射方案均選用RRTMG方案，陸面過程選用NOAH方案，邊界層方案選用YSU方案，積云對流參數(shù)化方案選用Kain-Fritsch方案。同化所需的分析數(shù)據(jù)即為提供模式初始場的NCEP全球?qū)α鲗臃治鰯?shù)據(jù)[16]。

表2 WRF模擬嵌套網(wǎng)格參數(shù)設(shè)置Table 2 Configuration of nested domains for WRF simulation

1.2 同化方案設(shè)置

根據(jù)筆者所在團隊的經(jīng)驗建議，設(shè)定3個采取不同同化方案的算例。其中，算例1(CTS)為控制試驗，不采用任何同化方案；算例2(AN)則在3個嵌套域內(nèi)均使用了AN同化方案，同化所需的分析數(shù)據(jù)同為提供模式初始場的NCEP全球?qū)α鲗臃治鰯?shù)據(jù)[18]，同化周期為6 h；算例3(SAN)則同時采用AN(同化周期為6 h)和SAN(同化周期為3 h)同化方案，SAN方案所需的近地面觀測數(shù)據(jù)采用NCEP全球近地面觀測氣象資料[19]。算例2、3的松弛系數(shù)均取默認值0.000 3 s-1。

1.3 中尺度WRF模擬結(jié)果

圖2給出了具有代表性意義的時間點下，對應(yīng)d03區(qū)域，AN方案模擬10 m離地高度處對應(yīng)的風(fēng)場演化云圖結(jié)果。在模擬時間內(nèi)，近地風(fēng)場風(fēng)速大致在0～7 m/s范圍內(nèi)，其中最大風(fēng)速出現(xiàn)在9月13日14時左右，結(jié)合當?shù)氐孛蔡卣骺梢?，高風(fēng)速區(qū)域主要集中在d03區(qū)域西部的海面部分，此范圍內(nèi)風(fēng)速可達5 m/s以上，隨著氣流深入內(nèi)陸，風(fēng)速逐漸減小，東部內(nèi)陸區(qū)域近地風(fēng)速基本在3 m/s以內(nèi)。從風(fēng)向來看，是以海向風(fēng)為主，氣流基本從西南海域方向去往東北內(nèi)陸方向。

將d03區(qū)域模擬風(fēng)場數(shù)據(jù)(水平方向共2 601個數(shù)據(jù)點)與NCEP氣象數(shù)據(jù)[16]進行對比，圖3為各時刻3個算例的10 m離地高度風(fēng)速模擬誤差云圖。由CTS模擬結(jié)果可見，正值誤差主要分布在d03西部沿海地區(qū)及東北部內(nèi)陸區(qū)域，并且此范圍內(nèi)誤差大小一般在1 m/s以上，而中部區(qū)域誤差主要以負值為主且誤差較小，整個模擬時域內(nèi)氣象塔鄰近范圍誤差不超過1 m/s；13日14時海面風(fēng)速誤差最為明顯，西部邊緣處誤差達到2 m/s。對比加載同化模塊(AN和SAN)后的誤差結(jié)果可見，F(xiàn)DDA技術(shù)對模擬結(jié)果進行了不同程度的糾偏，其中AN方案同化效果最為明顯，如13日2時，AN方案在東北角區(qū)的模擬誤差出現(xiàn)了大面積的降低，默認方案在此區(qū)域風(fēng)速偏差均在1 m/s以上，北部邊界處誤差甚至接近3 m/s，而AN方案將此范圍內(nèi)大部分區(qū)域的風(fēng)速誤差控制在1 m/s以下，北部少量塊狀區(qū)域最大誤差不超2 m/s。

選取深圳氣象梯度塔實測數(shù)據(jù)為參照，表3給出了3種數(shù)值模擬方法的模擬誤差比較。由表可見，WRF模擬結(jié)果與實測數(shù)據(jù)風(fēng)速值的均方根誤差(RMSE)大部分在1～2 m/s之間，其中AN方案誤差值低于其他兩個方案，相關(guān)系數(shù)則基本在0.7以上，表明數(shù)值模擬結(jié)果尤其是AN方案，與實測數(shù)據(jù)有較強的相關(guān)性；WRF模式在高空風(fēng)場的模擬均方根誤差小于近地面，表明WRF模式對高空處風(fēng)場的風(fēng)速模擬效果優(yōu)于近地面，其原因可能在于WRF模式采用的地形高程數(shù)據(jù)精度不足，且近地風(fēng)速的模擬易受復(fù)雜地形影響。

2 CFD降尺度風(fēng)場模擬

2.1 入口風(fēng)剖面公式擬合

選用前文2019年9月13日14時AN方案模擬的風(fēng)場結(jié)果作為降尺度CFD數(shù)值風(fēng)洞模型的入流邊界條件，根據(jù)模擬結(jié)果已知該時段氣象塔位置風(fēng)向角約為222°，因此選定圖4(a)所示的矩形區(qū)域ABCD作為CFD模擬區(qū)域(平面實尺尺寸為1 000 m×500 m)，中心圓點表示氣象塔所在位置，線段AB表示風(fēng)速入口，沿AB等間距布置11個點位作為入口風(fēng)速剖面提取點，具體點位布置如圖4(b)所示。

(a)2019年9月13日2時

(b)2019年9月13日8時

(c)2019年9月13日14時

(a)CTS，13日2時

(b)CTS，13日8時

(c)CTS，13日14時

(d)AN，13日2時

(e)AN，13日8時

(f)AN，13日14時

(g)SAN，13日2時

(h)SAN，13日8時

(i)SAN，13日14時

圖4(c)和4(d)分別給出了入口處11個提取點的506 m以下平均風(fēng)向角剖面及平均風(fēng)速剖面圖。由風(fēng)向角結(jié)果可見，11個點位的風(fēng)向角大小及其隨高度的變化規(guī)律基本一致，各點的350 m以下風(fēng)向角平均值約為222°；由圖4(d)可見，各點風(fēng)速剖面基本一致，表明了AB面對應(yīng)風(fēng)速具有較好的空間均勻性。如圖4(d)所示，模擬得到的平均風(fēng)剖面和規(guī)范規(guī)定的良態(tài)風(fēng)指數(shù)率剖面明顯不同，這可能與丘陵地形影響以及當時的風(fēng)氣候條件有關(guān)，后文將給出氣象塔實測結(jié)果的比較。

表3 各豎向?qū)语L(fēng)速模擬值均方根誤差及相關(guān)系數(shù)Table 3 RMSE and correlation coefficient of wind speed simulation values in each vertical layer

現(xiàn)取11個點位置處的空間平均值作為擬合入口風(fēng)剖面所需的數(shù)據(jù)點，擬合曲線如圖4(e)所示，對應(yīng)擬合公式見式(1)：

U=2.238 35+0.055 45z-4.454 22×10-4z2+

1.581 07×10-6z3-2.661 08×10-9z4+

1.701 51×10-12z5

(1)

式中：U為風(fēng)速值，m/s；z為海拔高度，m。

圖4 CFD計算域平面圖及入口風(fēng)剖面圖

2.2 復(fù)雜地形建模及網(wǎng)格劃分

圖5給出了CFD復(fù)雜地形建模及網(wǎng)格劃分示意圖。首先基于Arcgis平臺將ASTER GDEM地形高程數(shù)據(jù)轉(zhuǎn)化為三維空間坐標，并導(dǎo)入CATIA軟件進行實尺度山體建模(見圖5(a))，繼而交由ICEM軟件進行網(wǎng)格劃分，最終實尺計算域網(wǎng)格模型如圖5(b)所示，其長寬高為1000 m×500 m×500 m，坐標系原點為氣象塔所在位置，X=-500 m處為速度入口。算例采用非結(jié)構(gòu)化網(wǎng)格進行網(wǎng)格劃分，最大網(wǎng)格尺寸為20 m，近地面采用加密網(wǎng)格，如圖5(c)所示，對應(yīng)尺寸為2 m；在氣象塔位置處須進行網(wǎng)格加密以提高精度，故在此區(qū)域施加一圓柱形網(wǎng)格密度盒(見圖5(d))，其中心位置網(wǎng)格尺寸為1 m，并以1.2倍的延伸比向外擴散；網(wǎng)格總體數(shù)量為705萬。

(a)復(fù)雜地形建模

(b)整體計算流域及網(wǎng)格劃分

(c)地表網(wǎng)格

(d)流域中心網(wǎng)格密度盒

2.3 平衡態(tài)大氣邊界層驗證

平衡態(tài)大氣邊界層要求流域中湍流風(fēng)的速度剖面順風(fēng)向梯度為零[20]，參考Yang等[21]的建議對空域(即未加載地形模型的計算域)進行邊界層自保持性驗證。算例采用RANS方法進行求解，其中入口邊界采用式(1)生成入口平均風(fēng)速，各項參數(shù)設(shè)置見表4。

表4 CFD求解參數(shù)設(shè)置Table 4 Setting of CFD solution

圖6為空域Y=0 m截面自保持性驗證結(jié)果，由圖6(a)可見，風(fēng)速在整個流域沿順風(fēng)向基本保持恒定，圖6(b)顯示入口處(X=-500 m)、氣象塔位置(X=0 m)、出口處(X=500 m)對應(yīng)的風(fēng)廓線基本一致，由此判定數(shù)值模擬的邊界層風(fēng)速剖面自保持性良好，可滿足精細化模擬要求。

2.4 CFD模擬結(jié)果分析

將地形模型載入計算域，數(shù)值模型參數(shù)沿用表4 進行設(shè)置，圖7(a)給出了WRF與CFD降尺度耦合方法(下文簡稱CFD方法)模擬所得的氣象塔風(fēng)廓線，并將其與單一WRF模擬結(jié)果(WRF)、實測數(shù)據(jù)進行比較。

由圖7(a)可見：整體上，WRF與CFD方法模擬得到的風(fēng)剖面整體規(guī)律和氣象塔吻合，均異于規(guī)范中的良態(tài)風(fēng)指數(shù)率規(guī)律；對于50 m以下近地區(qū)域內(nèi)，單一WRF、CFD方法的模擬結(jié)果均略小于實測數(shù)據(jù)，而相對于CFD方法，WRF的結(jié)果反而與實測數(shù)據(jù)更為接近，這其中的誤差可能源自CFD平面風(fēng)場計算結(jié)果如圖7(b)-7(e)所示，分析可知：順風(fēng)向及橫風(fēng)向風(fēng)速在海拔150 m以上基本不發(fā)生變化，而在100 m以下受復(fù)雜地形影響，氣流在地表處因撞擊山體造成大量湍流渦旋運動，近地面風(fēng)速變化梯度明顯，山坡背風(fēng)面風(fēng)速明顯減弱，山坡頂處風(fēng)速沿豎向出現(xiàn)了明顯的加速特征；在100 m以下，風(fēng)速基本隨著離地高度的增加而增大，在100～150 m區(qū)間內(nèi)，風(fēng)速達到4.5 m/s以上，而隨著高度的繼續(xù)增加，風(fēng)速緩慢減小且趨于平緩，在450 m以上時，風(fēng)速降至3.5 m/s以下；對于水平截面的風(fēng)場分布，受風(fēng)場入口處山體隆起的影響，入口近地風(fēng)速掠過山頂后在背風(fēng)面出現(xiàn)了明顯的減速現(xiàn)象，導(dǎo)致兩側(cè)邊界鄰近范圍風(fēng)速大于中心風(fēng)速，而Y=250 m邊界處近地風(fēng)速略大于Y=-250 m 邊界處的近地風(fēng)速。

(a)Y=0 m截面風(fēng)速分布

(b)Y=0 m處入口、氣象塔及出口風(fēng)廓線比較

方法模擬區(qū)域的局限性，在接近流場壁面的區(qū)域模擬精度降低，且受入口處地形影響導(dǎo)致風(fēng)速下降；對于100 m以上高空風(fēng)場，WRF與CFD方法結(jié)果對應(yīng)的最大相對誤差為8.34%和5.07%，平均相對誤差為2.86%和1.46%，均方根誤差分別為0.20 m/s、0.13 m/s，其對應(yīng)相關(guān)系數(shù)分別為0.922 1、0.924 4，由此可見，WRF與CFD耦合技術(shù)使得高空風(fēng)場的數(shù)值模擬精度得到了進一步提升。

(a)氣象塔址處風(fēng)廓線模擬結(jié)果和實測結(jié)果比較

(b)X=-350 m，X=0 m，X=350 m截面風(fēng)速

(c)Y=0 m截面風(fēng)速分布

(d)海拔50 m處風(fēng)場

(e)海拔70 m處風(fēng)場

3 結(jié)論

本文選取深圳氣象塔周邊復(fù)雜地形區(qū)域，基于WRF與CFD進行了對該區(qū)域風(fēng)場的多尺度耦合模擬，并將模擬結(jié)果與氣象實測數(shù)據(jù)進行對比，得出如下結(jié)論：

(1)基于中尺度WRF模式，可重現(xiàn)復(fù)雜地形風(fēng)場的中尺度發(fā)展演化過程及整體風(fēng)場特性，其中AN同化方案模擬所得中尺度風(fēng)場結(jié)果與實測數(shù)據(jù)更為接近，且風(fēng)速模擬結(jié)果較CTS方案有很大提升，說明AN同化方案對WRF風(fēng)速模擬值有改進效果；AN方案各高度層對應(yīng)風(fēng)速模擬結(jié)果的最大均方根誤差不超1.53 m/s，最大相對誤差為8.12%，平均相對誤差為3.6%，表明該方法對中尺度風(fēng)場的模擬具有一定的準確性。

(2)從中尺度WRF模擬邊界層風(fēng)場數(shù)據(jù)中提取CFD流場仿真所需的入口邊界條件，并進行了平衡態(tài)大氣邊界層的驗證，運用RANS方法對氣象塔所在復(fù)雜地形區(qū)域?qū)崿F(xiàn)了精細化的降尺度數(shù)值模擬。模擬結(jié)果顯示，相對于僅采用WRF模擬，WRF-CFD多尺度耦合方法可較準確描述近地面風(fēng)場特性，如對于100 m以上高空風(fēng)場，其平均相對誤差為1.46%，風(fēng)速均方根誤差僅為0.13 m/s，對應(yīng)相關(guān)系數(shù)高達0.924 4，WRF-CFD多尺度耦合方法使數(shù)值模擬精度得到了進一步的提升。

(3)數(shù)值模擬和氣象塔實測結(jié)果顯示，對于本文研究的這一山地風(fēng)場，其平均風(fēng)剖面在選定的模擬時刻異于規(guī)范的良態(tài)風(fēng)指數(shù)率模型，也說明了在復(fù)雜地貌下邊界層風(fēng)場的復(fù)雜性，CFD風(fēng)場模擬結(jié)果部分解釋了這一現(xiàn)象。

未來將借助更先進的激光測風(fēng)雷達開展自然風(fēng)場的實測研究，結(jié)合本文基于多尺度模擬技術(shù)的湍流風(fēng)場精細化模擬技術(shù)和研究框架，擬為今后建筑結(jié)構(gòu)抗風(fēng)的精細化評估提供參考和奠定基礎(chǔ)。