邵應泉 , 張立杰 , 饒華炎

1. 珠海市國家氣候觀象臺， 廣東 珠海 519000 2. 深圳市國家氣候觀象臺， 深圳市氣象局， 廣東 深圳 518040

0 引 言

作為海洋資源的重要組成部分，海島與陸地自然資源一樣，是人類文明發(fā)展的重要支撐之一，也是人類文明向海洋延伸的踏板。在我國南海和東海區(qū)域，海島還是重要的國防屏障，對于保衛(wèi)國土完整和國家安全有重要意義。然而，海島通常也是生態(tài)系統(tǒng)極脆弱的區(qū)域，極易受到氣候變化和人類活動的影響(孫立廣等，2005；張浪等，2011)。因此，加深對海島的氣候特征的理解，對于保障海島的生態(tài)環(huán)境，利用海島的氣候資源有著重要意義。這其中，對于海島的風場特征的認識又特別重要，它不僅是充分認識海島資源稟賦、環(huán)境氣象條件與氣候背景的重要基礎，也是戰(zhàn)時指揮決策的重要科學依據。

然而，當前我國對于海島地區(qū)風場規(guī)律和路面過程特征的研究仍有不足，相關的科學認識也不甚明晰。特別的，相較于大陸地區(qū)已有的水—陸交界非均勻邊界層及陸面過程研究(吳楠等,2019；饒志娟等，2021)，在海島地區(qū)開展的研究則相對更少，且組織開展觀測實驗難度也更大。盡管有一些學者在海島地區(qū)組織了地面風場觀測實驗(黃浩輝等，2007)，但由于只能在有限的點位上部署儀器開展觀測，所獲得的研究成果對于整個海島的代表性也存在一定不足，而關于風資源利用、環(huán)境影響評價等工作，需要能覆蓋整個島嶼的風場數(shù)據。海島作為凸出于海面的陸地，地形通常較為復雜，其地面風場的非均勻性也會較為突出，在這種背景下，依靠數(shù)值模擬手段分析海島的風場特征特別重要。

目前依靠數(shù)值模擬手段分析海島的風場特征有著較大困難——在我國的領海，除了臺灣島和海南島外，大多數(shù)島嶼的空間尺度較小，多為1 km以內量級或更小的島嶼。而大氣科學中常用的數(shù)值模擬工具通常為中尺度數(shù)值模式，這些模式的空間網格距通常在1 km數(shù)量級，盡管也有利用數(shù)值模式展開分辨率為100 m的模擬研究，但總體上由于模式在數(shù)值求解時多采用有限差分法，因而對過于復雜的地形通常需要采用平滑處理，從而使得模式中的地形與真實地形存在較大差距(Li et al，2010)。

事實上，國外有不少研究在開展海島風場研究時，采用了計算流體力學(Computational Fluid Dynamics，CFD)工具。與中尺度模式相比，CFD通常是基于有限體積法進行數(shù)值求解，因此可適用于任意復雜的幾何形體，自然也適用于復雜地形(李磊等，2010)。英國科學家于20世紀80年代在Askervein海島地區(qū)開展過一系列觀測實驗和模擬實驗，以研究復雜地形邊界層的風特性(Taylor and Reunissen，2011)。歐盟科學家則針對丹麥的Bolund島開展過模擬研究，以驗證CFD模型中的一些模型(Berg et al，2011；Chaudhari et al，2016)。盡管這些研究已經取得的了一定進展，但關于海島風場模擬的研究總體上仍然較少，特別是在我國，幾乎未曾見到針對海島特別是小型海島風場的模擬研究。這是因為關于海島的高精度高程數(shù)據并不容易獲取，從而為構建模型CFD計算模型造成了障礙。

基于上述原因，文中提出一種技術思路，直接利用SRTM(Shuttle Radar Topography Mission)資料提取海島地形數(shù)據，建立精細化的CFD計算模型，并利用模型進行數(shù)值模擬試驗。該項技術有望為海島風資源評估、風電場選址、污染擴散等領域的研究提供參考。

1 資料與技術路線

SRTM數(shù)據主要是由美國航空航天局(NASA)和美國國防部國家測繪局(NIMA)聯(lián)合測量獲取的地形數(shù)據集(湯國安等，2010)，英文全稱為“Shuttle Radar Topography Mission”，可譯作航天飛機雷達地形測繪使命(下載地址：http://srtm.csi.cgiar.org)。在該數(shù)據的官方網站上，將全球劃分為若干個區(qū)域，每個區(qū)域為一個獨立的文件，以“TIF”格式存儲。文中使用的數(shù)據文件為“SRTM_59_08.tif”文件，該文件包含(109.9995838°—115.0004172°E,19.9995838°—25.0004171°N)范圍內的所有區(qū)域，空間分辨率為0.00083°，在赤道地區(qū)約為90 m。我國珠江口附近的島嶼均包含在該文件中，文中特別截取了臨近珠海東澳島的區(qū)域進行試驗。

使用的CFD軟件為FLUENT通用計算流體力學軟件，它的前處理器采用了GAMBIT，用于生成計算機輔助設計(CAD)模型并進一步劃分網格，形成可被FLUENT讀取的解域文件。具體技術路線如圖1所示：1) 下載獲取SRTM原始數(shù)據文件；2) 利用地理信息系統(tǒng)(GIS)工具讀取原始文件，通過交互式方式選取并分割海島地形高程數(shù)據，生成ASCII碼格式的高程數(shù)據文本文件；3) 用FORTRAN語言編寫GAMBIT腳本生成器，該生成器的功能主要是讀取高程數(shù)據文本文件，并生成腳本，腳本可被GAMBIT讀取自動執(zhí)行；4) 利用GAMBIT讀取生成的腳本文件并自動執(zhí)行，形成高程點陣每個點陣為笛卡爾坐標系下坐標為(x, y, z)的點，然后通過交互式操作生成點陣的包絡面，在建立包絡面的基礎上進一步生成體并劃分網格；5) 輸出.msh文件，該文件可被FLUENT等CFD軟件識別和讀取。

圖1 SRTM數(shù)據處理技術路線

在生成網格文件之后，利用FLUENT進行風場模擬，其基本方程如式(1)—(3)所示：

(1)

(2)

(3)

2 結果與分析

2.1 地形數(shù)據與數(shù)據提取

利用GIS軟件讀取“SRTM_59_08.tif”文件，其可視化界面如圖2a所示。由圖可見，該文件包含了幾乎整個廣東省及南海的部分區(qū)域，北至南嶺山地，南至雷州半島，珠江口附近的小島均被包含在文件中。通過逐級放大，可鎖定目標小島(如圖2b中黑色方框所示)，裁取包含小島的數(shù)據片段，即可獲得生成CAD模型及相關腳本文件所必須的地形高程數(shù)據文件。

圖2 鎖定模擬區(qū)域及提取所需數(shù)據(a.地形文件；b.東澳島)Fig. 2 Simulation area and the extracted data (a. Terrain file; b. Dong’ao island)

2.2 CAD模型與網格文件

利用FORTRAN程序編程得到GAMBIT腳本生成器。運行該生成器，可以根據上一步得到的數(shù)據文件生成一個腳本，腳本的每一行為一個自動操作命令，按照給定的坐標生成一個“點”。用于生成東澳島CAD模型的腳本共有5 440行，在GAMBIT中運行腳本，相當于批處理執(zhí)行生成“點”的操作，即建立了一個包含5 440個“點”的點陣，如圖3a所示。進一步基于點陣，利用GAMBIT的面生成功能生成一個包絡面，即可以較好地展現(xiàn)出海島的地形特征(圖3b)。

圖3 GAMBIT生成的點陣(a)和基于點陣建立的包絡面(b)Fig. 3 Dot matrix generated by GAMBIT (a) and envelope surface based on the matrix (b)

在包絡面基礎上，依次生成4個側面和1個頂部的平面，并組合成一個“體”，即可作為模擬區(qū)域。頂部的高度一般應設置為地形高度的10倍，考慮到試驗選取的東澳島高度不到200 m，故解域的頂高設置為2 000 m。對比東澳島CAD模型地面高程分布與衛(wèi)星圖片，可知文中建立得到的CAD模型無論從島的海岸輪廓還是從地形起伏的特征來看，都較為一致。由此可見，盡管SRTM數(shù)據是開源數(shù)據，但基于該數(shù)據建立的海島模型也很準確，足以為海島風場的模擬研究奠定基礎。

生成解域之后，須根據具體需求設置網格。文中使用了非均勻網格的設置，以圖4a中的十字線為中心，分別向東西和南北兩個方向逐漸拉伸網格距。東西和南北方向的網格數(shù)分別為200和160，水平方向上設置為30層，其中最底層的高度約為離地2 m左右。解域內下墊面的網格設置如圖5所示。

圖4 東澳島CAD模型的地面高程分布(a)和東澳島衛(wèi)星遙感圖像(b)Fig. 4 Ground elevation distribution of Dong'ao island by CAD model (a) and satellite remote sensing image of theisland (b)

圖5 東澳島風場模擬實驗的地面網格設置

采用非均勻網格，是希望在東澳島的東側獲得更高分辨率的數(shù)據，該側為氣流的來流方向，對于研究迎風側的氣流爬升等細節(jié)問題有幫助。而對于暫不關注的其他區(qū)域，采用較稀疏的網格設置，有利于節(jié)約計算資源，獲得更快的模擬速度。

2.3 風場模擬結果

東澳島東側設有一個自動氣象站，文中利用該站在2020年9月15日18:30左右監(jiān)測到數(shù)據作為輸入參數(shù)。該站在10 m高度處測得的2 min平均風速為5.7 m/s，入流風向角為69°(大致為東北偏東風向)。根據該觀測數(shù)據，設入流風速的垂直廓線如式(4)所示：

(4)

式中，U(z)為z高度上的風速，U10為10 m高度的風速，z為距地面高度，0.2為指數(shù)風廓線在中性層結下的指數(shù)取值。將風速分別在東西和南北方向上分解，得到u分量和v分量，再將2個分量在垂直方向上的變化用FLUENT自帶的用戶自定義函數(shù)(UDF)編寫出來。在正式模擬時，將事先編寫的UDF載入模型，并在入流邊界上賦值。

模擬分兩步進行，第一步為穩(wěn)態(tài)模擬，即不在模擬中加入時間偏導項。如前所述，在穩(wěn)態(tài)模擬時使用雷諾平均的湍流閉合方案，realizable k-ε方案。待穩(wěn)態(tài)模擬達到收斂以后，開啟時間偏導項，并將湍流閉合方案設置為大渦模擬(LES)方案，積分3 600 s后，得到最終的模擬結果。

圖6給出了模擬得到的距地面10 m高度的風矢量分布。由圖6a可知，模擬結果較好的反映了復雜地形對風場的影響，在海拔較高的地方風速均較大，氣流在爬越山峰的過程中，風速有明顯增大的趨勢，這與Li等(2014)之前在深圳排牙山地區(qū)的模擬結果類似。此外，圖6a還表明，非均勻設置的網格有效地捕捉到了所關心區(qū)域的風場細節(jié)，在關注的東澳島東側地區(qū)，氣流在接觸島嶼時受地形阻礙損失動能形成的相對小風區(qū)較為明顯。而在網格較為稀疏的其他地區(qū)，這一現(xiàn)象并不明顯。圖6b還給出了局部的模擬結果，由圖可觀察到更多的細節(jié)，特別是在背風側，除了風速顯著減小之外，風向也受地形影響出現(xiàn)了明顯變化，在島嶼西側的部分區(qū)域，風向已由東北偏東變?yōu)榱似憋L，這是氣流爬越和繞開島嶼后繼續(xù)前進的表現(xiàn)。

圖6 2020年9月15日18:30東澳島風矢量(填色，單位：m/s)水平分布全貌(a)和局部(b)模擬結果Fig. 6 Overall view (a) and local view (b) of horizontal distribution of wind vector (shaded；unit:m/s) on Dong’ao island at 18:30 BT on September 15, 2020

圖7進一步給出了模擬區(qū)域內沿主導風向垂直剖面上的風場分布，由圖可見氣流爬越島嶼的整個過程。在垂直方向上，風速隨高度增高而增加的趨勢十分明顯。而在島嶼的山體背后，由于氣流爬越山體損失部分動能，風速明顯比山前要小。綜上所述，利用文中提出的技術思路，可以較好地實現(xiàn)島嶼邊界層風場的模擬，且模擬結果在物理上可以得到合理解釋，這主要體現(xiàn)在風速的增強、衰減及風向的轉變、擾動都符合預期，與地形的起伏有很好的對應關系，并與已有的研究結果相吻合。

圖7 2020年9月15日18:30東澳島風矢量(填色，單位：m/s)垂直分布模擬結果

4 結論與展望

文中提出了利用SRTM資料構建海島的CFD計算模型的技術方法，并以珠江口的東澳島為對象，進行了模型構建試驗及邊界層風場的數(shù)值模擬試驗。實驗結果表明：

1) SRTM資料作為開源的資料，可有效解析空間尺度在1 km內數(shù)量級的海島地形；

2) 基于SRTM資料，可以構建用于CFD模擬的計算模型，并較為真實地刻畫出島嶼的地形；

3) 利用計算模型進行的模擬試驗表明，島嶼的邊界層風場受島嶼地形影響十分明顯，島嶼上空10 m高度的風場隨地形增加而明顯增加，在島嶼的背風側，風力明顯減弱，且風向有十分明顯的變化。

必須指出，本研究仍然只是初步研究，重心放在了計算模型的構建上，對于島嶼邊界層風場結構的科學分析比較粗淺。同時，由于目前觀測資料的限制，模擬數(shù)據與觀測數(shù)據的對比還未曾涉及，這些都是未來研究需要進一步彌補的地方。盡管如此，本文從技術上打通了從SRTM數(shù)據到CFD模擬的全流程，這些技術未來在海島風資源評估、風電場選址、污染擴散等領域可以發(fā)揮作用。