董浩天，陶 韜，杜曉慶

(1.上海大學(xué)，上海 200444；2.東京工業(yè)大學(xué)，日本 東京 152-8550)

0 引 言

我國東南沿海常年受到西北太平洋臺風(fēng)災(zāi)害的影響[1-2]。據(jù)國家應(yīng)急管理部統(tǒng)計，僅2018年前三季度全國臺風(fēng)等風(fēng)致災(zāi)害就造成4 712.2萬人次受災(zāi)，203人死亡，直接經(jīng)濟(jì)損失854.2億元。隨著計算機(jī)硬件和軟件的發(fā)展，近年來中尺度氣象模式常用于臺風(fēng)過程的模擬和預(yù)測。WRF（Weather Research andForecasting model）開源中尺度氣象模式采用了完全可壓縮非靜力平衡歐拉方程，并提供了大量物理參數(shù)化方案和數(shù)值方法選項，在大氣邊界層模擬和預(yù)報上明顯優(yōu)于傳統(tǒng)大尺度氣象模式[3]。Gu等[4]采用WRF模式模擬了2002年臺風(fēng)Rusa的中心路徑和風(fēng)速分布。Tse等[5]采用WRF模擬了臺風(fēng)過程近地風(fēng)場，發(fā)現(xiàn)WRF可得到較好的平均風(fēng)剖面和風(fēng)速時程結(jié)果，但近地風(fēng)速誤差顯著。Jiménez等[6]通過在控制方程中加入由地形復(fù)雜度定義的源項，提高了WRF模擬亞格子地形特征的能力。總的來說，WRF可以一定程度上預(yù)測/再現(xiàn)臺風(fēng)路徑，并較好地模擬近地10 min平均風(fēng)速時程；但由于控制方程形式、數(shù)值過濾、貼體網(wǎng)格和邊界條件等原因，其對臺風(fēng)過程復(fù)雜地形近地湍流特征的模擬效果仍難以滿足風(fēng)工程的研究需求[7]。

復(fù)雜地形風(fēng)場研究，特別是我國東南沿海地區(qū)復(fù)雜地形臺風(fēng)過程風(fēng)場的研究，對災(zāi)害預(yù)防和風(fēng)資源開發(fā)有著重要意義[8-10]。而同平坦地形相比，峽谷、高山、海峽等山區(qū)地形的特征尺度與大氣邊界層高度相近，易發(fā)生越山風(fēng)、峽谷風(fēng)和特征湍流等現(xiàn)象。隨著湍流模型、網(wǎng)格技術(shù)、數(shù)值算法和計算機(jī)硬件的發(fā)展，計算流體動力學(xué)（CFD）方法已廣泛應(yīng)用于風(fēng)工程領(lǐng)域并部分替代風(fēng)洞試驗，新版抗風(fēng)規(guī)范[11]已對采用CFD虛擬風(fēng)洞試驗方法進(jìn)行橋梁抗風(fēng)研究給出了指南。Tamura[12]指出，大渦模擬方法（Large-eddy Simulation，LES）在非定常流動細(xì)節(jié)和計算資源消耗之間可達(dá)到較好的平衡，適用于結(jié)構(gòu)抗風(fēng)設(shè)計、復(fù)雜地形風(fēng)場、城市風(fēng)環(huán)境等風(fēng)工程課題。Cao等[13]采用LES和浸入邊界法模擬了湍流邊界層風(fēng)流經(jīng)粗糙二維小山。復(fù)雜地形風(fēng)場具備多尺度特征，即存在不同尺度、不同頻率成分的大氣流動。抗風(fēng)規(guī)范[11]采用準(zhǔn)定常假設(shè)，以10 min時距風(fēng)速測量和規(guī)范風(fēng)剖面、風(fēng)譜構(gòu)建平穩(wěn)來流條件，可以較好描述小尺度流動現(xiàn)象，但忽視了橋址風(fēng)場的多尺度特征。

為了克服中尺度模式在復(fù)雜地形近地湍流特征模擬方面的不足，Perivolaris等[14]學(xué)者提出了中尺度氣象模式同小尺度CFD結(jié)合的多尺度數(shù)值模擬方法。如圖1[15]所示，通過WRF嵌套網(wǎng)格求解大、中尺度三維風(fēng)場低頻信息；通過人工湍流生成技術(shù)補充高頻風(fēng)速成分，生成多尺度湍流風(fēng)、溫場，并作為邊界條件輸入CFD計算，求解近地風(fēng)速時程。高頻風(fēng)速成分的再生成一直是多尺度數(shù)值模擬的難點。Wyszogrodzki等[16]直接將WRF結(jié)果經(jīng)時間空間插值后輸入LES城市地貌計算域，其上游地表起到了“數(shù)值粗糙元”的作用，湍流風(fēng)剖面沿流向自然發(fā)展，但這一方法要求很長的上游計算域長度。Nakayama等[17]通過循環(huán)入流技術(shù)在WRF平均風(fēng)剖面基礎(chǔ)上進(jìn)行人工湍流生成，降低了計算資源消耗。

圖1 多尺度風(fēng)場模擬技術(shù)路線Fig.1 A sketch of the multi-scale simulation of a wind field over a complex terrain

本文以2012年臺風(fēng)啟德登陸廣東東海島為背景，采用多尺度風(fēng)場模擬方法，研究了臺風(fēng)路徑、近地風(fēng)速時程和風(fēng)剖面等，并同實測結(jié)果驗證。研究內(nèi)容可供東南沿海臺風(fēng)影響下的臺風(fēng)災(zāi)害預(yù)報、結(jié)構(gòu)抗風(fēng)設(shè)計和風(fēng)力發(fā)電選址等參考。

1 臺風(fēng)模擬資料

1.1 臺風(fēng)啟德概況

1213號臺風(fēng)啟德（Kai-tak）于2012年8月13日在菲律賓以東洋面加強(qiáng)成為第13號熱帶風(fēng)暴；8月15日21時加強(qiáng)為臺風(fēng)；8月17日4時30分前后在廣東省湛江市登陸，登陸時中心附近最大風(fēng)力達(dá)到13級；之后越過雷州半島，進(jìn)入北部灣并在越南再次登陸，最后于8月18日逐漸消散。圖2顯示了啟德臺風(fēng)的中心路徑。

圖2 臺風(fēng)啟德中心路徑（來源：日本氣象廳網(wǎng)站）Fig.2 The central route of typhoon Kai-tak 1213

1.2 氣象觀測資料

測站觀測資料主要用于對近地風(fēng)場模擬結(jié)果進(jìn)行驗證。廣州省湛江市東海島布置的測風(fēng)塔記錄了該次臺風(fēng)過境時的風(fēng)速相關(guān)信息。東海島緊鄰雷州半島，測風(fēng)塔的位置為（21° 0' 44.90' N，110° 31' 33.88' E），該位置位于臺風(fēng)路徑的北測。如圖3所示，測風(fēng)塔儀器安置高度分別在距離地表100 m、70 m、50 m、30 m高度處。測風(fēng)設(shè)備為NRG #40風(fēng)速計和NRG #200p風(fēng)向傳感器，輸出數(shù)據(jù)包括風(fēng)速和風(fēng)向的10 min平均值、標(biāo)準(zhǔn)差、最大值和最小值。

圖3 東海島測風(fēng)塔Fig.3 The wind observation tower at Donghai Island

1.3 氣象再分析資料

氣象再分析資料為多尺度數(shù)值模擬提供初值條件和邊界條件更新。本次對臺風(fēng)啟德模擬采用的再分析資料為美國國家環(huán)境預(yù)報中心（NECP）提供的FNL全球分析資料（Final Operational Global Analysis）。數(shù)據(jù)解析時間每6 h更新一次，格點分辨率為1°，包含了地表26個標(biāo)準(zhǔn)等力層（100～1 kPa）、地表邊界層、對流層頂?shù)囊匦畔ⅰ?/p>

1.4 地形地貌資料

在地表數(shù)據(jù)方面，對于WRF計算的外層d01和d02區(qū)域采用了美國地質(zhì)勘探局（USGS）提供的1'高精度高程數(shù)據(jù)和地表覆蓋數(shù)據(jù)；對于WRF內(nèi)層的d03區(qū)域和LES計算區(qū)域，則采用國家基礎(chǔ)地理信息中心提供的30 m精度ASTER GDEM地形數(shù)據(jù)和GLC地表覆蓋數(shù)據(jù)。

2 多尺度模擬方法

2.1 WRF計算參數(shù)

中尺度數(shù)值計算采用Advanced Research WRF（ARW，即WRF用于科學(xué)研究的計算內(nèi)核[3]）。ARW為基于全可壓縮非大氣靜力平衡方程，采用貼地的壓力定義的高度坐標(biāo)系，模式頂?shù)膲毫? kPa。WRF的計算域和網(wǎng)格參數(shù)如圖4和表1所示。豎向共61層網(wǎng)格，其中在觀測點位置最底部的10層高度分別為11.23、13.82、15.56、18.17、20.80、22.57、26.08、27.88、29.70、33.28 m。水平方向采用三層雙向嵌套網(wǎng)格，即d01、d02和d03，最外層d01包括了臺風(fēng)啟德登陸過程的區(qū)域，最內(nèi)側(cè)d03則主要覆蓋了東海島和周邊地形。每層網(wǎng)格的格點數(shù)、間距和時間步長如表1所示，時間積分采用了Runge-Kutta方法。物理參數(shù)化方案方面，采用了WSM6微物理方案、Kain-Fritsch卷積云方案、Monin-Obukhov表面層方案、YSU大氣邊界層方案、RRTMG輻射方案和Unified Noah陸面層方案。

圖4 WRF三層嵌套區(qū)域和臺風(fēng)啟德中心路徑結(jié)果Fig.4 Nested domains and typhoon routes

表1 啟德模擬三層嵌套網(wǎng)格參數(shù)設(shè)置Table 1 Parameters in the WRF simulation

2.2 循環(huán)入流人工湍流生成技術(shù)

采用循環(huán)入流技術(shù)[18-19]進(jìn)行高頻風(fēng)速成分的補充。如式（1）所示，空邊界層數(shù)值風(fēng)洞的入口速度分布(uin,vin,win)由WRF計算得到的平均速度剖面u0(z)和入口下游1倍邊界層高度處循環(huán)斷面速度(ure,vre,wre)的脈動部分組成，其中變量上劃線代表時間平均值。

2.3 LES計算參數(shù)

小尺度LES風(fēng)場計算的網(wǎng)格如圖5所示。計算域的頂部采用對稱邊界，海拔高度統(tǒng)一為2 143 m。下部邊界為無滑移壁面，海拔高度從0～78.28 m之間變化。北側(cè)和南側(cè)邊界均采用對稱邊界條件，計算域南北方向?qū)?00 m。東側(cè)采用2.2節(jié)得到的速度時程作為速度入口條件，西側(cè)則為出流條件，計算域東西方向長5 000 m。計算中主要模擬了來流臺風(fēng)邊界層從東側(cè)海平面向西越過沿海丘陵地區(qū)的過程，來流的平均方向為正西。水平方向網(wǎng)格劃分，東西向Δx和南北向 Δz均為 Δx= Δz= 10 m，共 500 × 80個網(wǎng)格。豎向共劃分150個網(wǎng)格，其中第一地表第一層網(wǎng)格高度Δymin= 0.143 m。六面體網(wǎng)格數(shù)量共計600萬。

圖5 小尺度LES計算網(wǎng)格示意圖Fig.5 The computational domain and boundary conditions for LES

LES部分入口的平均速度采用WRF內(nèi)層d03網(wǎng)格計算得到UTC時間2012年8月17日8：30附近的10 min平均速度剖面；風(fēng)速脈動部分則通過2.3節(jié)的循環(huán)入流方法生成，入口風(fēng)速剖面如圖9和圖10所示。計算中采用準(zhǔn)定常假設(shè)，即入口平均風(fēng)速分布不隨時間變化。邊界層頂?shù)臅r均風(fēng)速固定為Uref=27.82 m/s，邊界層高度固定約為δ99= 600 m。在計算中通過改變空氣黏性，以邊界層高度計的雷諾數(shù)約為5.7 × 104。LES計算的時間步長固定為0.256 8 s，計算中克朗數(shù)保持在0.4以下。采樣時間從t=1 027 s開始，此時以Uref= 27.82 m/s計的流體流經(jīng)流場5.7遍；采樣時長為600 s，對應(yīng)流經(jīng)遍數(shù)為3.3遍。

3 結(jié)果與討論

3.1 WRF臺風(fēng)過程模擬結(jié)果

圖4給出了WRF模擬得到的臺風(fēng)中心路徑時程結(jié)果和日本氣象廳提供的臺風(fēng)中心路徑觀測結(jié)果的對比。兩者的時距均為6 h，從UTC時間2012年8月16日20：00至17日20：00，方向從東南到西北。受到初值偏差的影響[20]，d01初始段臺風(fēng)路徑誤差較大，達(dá)到60 km左右。經(jīng)過網(wǎng)格雙向嵌套以后，內(nèi)層d03網(wǎng)格覆蓋區(qū)域的臺風(fēng)路徑精度較高，且明顯優(yōu)于外層d01和d02網(wǎng)格區(qū)域。

圖6和圖7給出了東海島測站位置距離地表100 m高度處水平風(fēng)速風(fēng)向24 h時程的WRF結(jié)果與觀測值對比。其中WRF和觀測結(jié)果均為10 min平均值。時間按照UTC標(biāo)準(zhǔn)時間給出。WRF和觀測得到的最大風(fēng)速分別為30.3 m/s和33.6 m/s，數(shù)值較為接近，且均出現(xiàn)在17日12:30時附近。從17日11:00-12:00時范圍的風(fēng)速低谷區(qū)域以及對應(yīng)的風(fēng)向大致180度的變化來看，WRF可以較好地模擬臺風(fēng)眼過境影響下的風(fēng)速、風(fēng)向變化，但一定程度上高估了臺風(fēng)眼過境的時長?？偟膩碚f，WRF可以較好地模擬臺風(fēng)登陸過程的大氣邊界層風(fēng)速時程，但模擬精度仍有待提高。

圖6 東海島測站100 m高度10 min平均風(fēng)速Fig.6 10-min averaged wind speed at 100 m above the ground

圖7 東海島測站100 m高度10 min平均風(fēng)向Fig.7 10-min averaged wind direction at 100 m above the ground

圖8中藍(lán)色點為LES計算得到的測站位置（x=0 m，z= 0 m）處地表10 m高度處x向風(fēng)速和y向風(fēng)速時程，而紅色線為通過WRF平均風(fēng)速+循環(huán)入流方法得到的入口（x= 2 500 m，z= 0 m）10 m高度處對應(yīng)的風(fēng)速時程。表2進(jìn)一步給出了統(tǒng)計得到的10 min平均風(fēng)速和湍流強(qiáng)度，其中湍流強(qiáng)度Iu=σu/Uref,Iv=σv/Uref，σu和σv分別為流向和豎向速度的標(biāo)準(zhǔn)差。可見，受到上游地形的影響，對于10 m高度處流向速度，測站處較入口平均速度下降8%，湍流強(qiáng)度則增加1倍左右；對于豎向速度，測站處出現(xiàn)0.73 m/s向下的平均速度，湍流強(qiáng)度稍弱于入口位置。圖9給出了入口和測站10 m高度東西向風(fēng)速的功率譜密度，可見測站位置的風(fēng)速高頻成分明顯多于入口。

圖8 測站位置和入口處10 m高度東西向風(fēng)速u和南北向風(fēng)速v時程Fig.8 Wind speed histories at 10 m above the ground

表2 測站位置和入口處10 m高度風(fēng)速統(tǒng)計值Table 2 Wind speeds at 10 m above the ground at the inlet and observation tower

圖9 測站位置和入口處10 m高度東西向風(fēng)速功率譜密度Fig.9 The power spectra density of u at 10 m above the ground

3.2 LES復(fù)雜地形風(fēng)場模擬結(jié)果

圖10和圖11給出了平均流向速度um和流向速度標(biāo)準(zhǔn)差σu的垂直剖面分布從上游到下游的變化。從地形來看，測站位置（x= 0 m）處于山頂偏背風(fēng)側(cè)。從入口（x= 2 500 m）開始，大約1 300 m左右的距離為海面。此外x= 1 000 m和x= -1 000 m兩個位置則分別位于測站所處的山頂?shù)纳嫌魏拖掠?。從平均風(fēng)剖面來看，測站處存在一定的越山加速現(xiàn)象，也即近地平均風(fēng)速的提高；LES的結(jié)果比WRF得到的近地風(fēng)速稍高，兩者均高于觀測值。從脈動風(fēng)剖面來看，受到地表流動分離的影響，LES得到的近地脈動風(fēng)速遠(yuǎn)高于WRF的結(jié)果，更接近觀測值；但50 m以上的脈動風(fēng)速下降明顯。此外，觀察到入口的脈動風(fēng)速剖面在200～300 m段以及500～600 m段存在較大的值。導(dǎo)致LES模擬誤差的主要原因可能有以下三點：1）WRF計算得到的平均風(fēng)剖面誤差傳遞到LES計算中；2）LES復(fù)雜地形風(fēng)場計算和循環(huán)入流湍流生成段采用了相同且較低的雷諾數(shù)，與實際大氣差異較大；3）公式（1）未對湍流向邊界層頂?shù)陌l(fā)展進(jìn)行阻尼約束。

圖10 上游（右）到下游（左）平均風(fēng)速剖面變化（z = 0 m位置）Fig.10 Profiles of the mean wind speeds (z = 0 m)

圖11 上游（右）到下游（左）脈動風(fēng)速剖面變化（z = 0 m位置）Fig.11 Profiles of the fluctuating wind speeds (z = 0 m)

4 結(jié) 論

采用WRF中尺度氣象模式和小尺度大渦模擬結(jié)合的復(fù)雜地形風(fēng)場多尺度模擬方法，研究了2012年臺風(fēng)啟德在廣東東海島登錄過程的臺風(fēng)路徑、100 m高度10 min平均風(fēng)速風(fēng)向時程、10 m高度近地風(fēng)速時程、平均和脈動風(fēng)剖面等結(jié)果，主要結(jié)論如下：

1）多尺度風(fēng)場模擬方法可實現(xiàn)臺風(fēng)影響下沿海地區(qū)復(fù)雜地形風(fēng)場時程的模擬，該方法不依賴規(guī)范風(fēng)剖面假定或現(xiàn)場實測，對沿海地區(qū)臺風(fēng)災(zāi)害的預(yù)報和評估具有較大的現(xiàn)實意義。

2）臺風(fēng)登陸過程中受到沿海復(fù)雜地形的作用，近地順流向湍流強(qiáng)度增加1倍左右，高頻風(fēng)速成分顯著增加，同時出現(xiàn)了0.73 m/s向下的平均豎向風(fēng)速。因此在臺風(fēng)登陸過程中，復(fù)雜地形局部可能產(chǎn)生遠(yuǎn)強(qiáng)于海面上的近地陣風(fēng)，對抗風(fēng)安全性產(chǎn)生不利影響。

3）后續(xù)研究可以從改善WRF風(fēng)剖面模擬精度、提高雷諾數(shù)和修改循環(huán)湍流生成方法等三個方面入手，進(jìn)一步提高近地風(fēng)場的模擬精度。

致謝：感謝中國氣象局、JMA、NECP、USGS、國家基礎(chǔ)地理信息中心等機(jī)構(gòu)提供的數(shù)據(jù)。