基于SRTM资料的小海岛CFD计算模型构建——以东澳岛为例

邵应泉 , 张立杰 , 饶华炎

1. 珠海市国家气候观象台， 广东 珠海 519000 2. 深圳市国家气候观象台， 深圳市气象局， 广东 深圳 518040

0 引 言

作为海洋资源的重要组成部分，海岛与陆地自然资源一样，是人类文明发展的重要支撑之一，也是人类文明向海洋延伸的踏板。在我国南海和东海区域，海岛还是重要的国防屏障，对于保卫国土完整和国家安全有重要意义。然而，海岛通常也是生态系统极脆弱的区域，极易受到气候变化和人类活动的影响(孙立广等，2005；张浪等，2011)。因此，加深对海岛的气候特征的理解，对于保障海岛的生态环境，利用海岛的气候资源有着重要意义。这其中，对于海岛的风场特征的认识又特别重要，它不仅是充分认识海岛资源禀赋、环境气象条件与气候背景的重要基础，也是战时指挥决策的重要科学依据。

然而，当前我国对于海岛地区风场规律和路面过程特征的研究仍有不足，相关的科学认识也不甚明晰。特别的，相较于大陆地区已有的水—陆交界非均匀边界层及陆面过程研究(吴楠等,2019；饶志娟等，2021)，在海岛地区开展的研究则相对更少，且组织开展观测实验难度也更大。尽管有一些学者在海岛地区组织了地面风场观测实验(黄浩辉等，2007)，但由于只能在有限的点位上部署仪器开展观测，所获得的研究成果对于整个海岛的代表性也存在一定不足，而关于风资源利用、环境影响评价等工作，需要能覆盖整个岛屿的风场数据。海岛作为凸出于海面的陆地，地形通常较为复杂，其地面风场的非均匀性也会较为突出，在这种背景下，依靠数值模拟手段分析海岛的风场特征特别重要。

目前依靠数值模拟手段分析海岛的风场特征有着较大困难——在我国的领海，除了台湾岛和海南岛外，大多数岛屿的空间尺度较小，多为1 km以内量级或更小的岛屿。而大气科学中常用的数值模拟工具通常为中尺度数值模式，这些模式的空间网格距通常在1 km数量级，尽管也有利用数值模式展开分辨率为100 m的模拟研究，但总体上由于模式在数值求解时多采用有限差分法，因而对过于复杂的地形通常需要采用平滑处理，从而使得模式中的地形与真实地形存在较大差距(Li et al，2010)。

事实上，国外有不少研究在开展海岛风场研究时，采用了计算流体力学(Computational Fluid Dynamics，CFD)工具。与中尺度模式相比，CFD通常是基于有限体积法进行数值求解，因此可适用于任意复杂的几何形体，自然也适用于复杂地形(李磊等，2010)。英国科学家于20世纪80年代在Askervein海岛地区开展过一系列观测实验和模拟实验，以研究复杂地形边界层的风特性(Taylor and Reunissen，2011)。欧盟科学家则针对丹麦的Bolund岛开展过模拟研究，以验证CFD模型中的一些模型(Berg et al，2011；Chaudhari et al，2016)。尽管这些研究已经取得的了一定进展，但关于海岛风场模拟的研究总体上仍然较少，特别是在我国，几乎未曾见到针对海岛特别是小型海岛风场的模拟研究。这是因为关于海岛的高精度高程数据并不容易获取，从而为构建模型CFD计算模型造成了障碍。

基于上述原因，文中提出一种技术思路，直接利用SRTM(Shuttle Radar Topography Mission)资料提取海岛地形数据，建立精细化的CFD计算模型，并利用模型进行数值模拟试验。该项技术有望为海岛风资源评估、风电场选址、污染扩散等领域的研究提供参考。

1 资料与技术路线

SRTM数据主要是由美国航空航天局(NASA)和美国国防部国家测绘局(NIMA)联合测量获取的地形数据集(汤国安等，2010)，英文全称为“Shuttle Radar Topography Mission”，可译作航天飞机雷达地形测绘使命(下载地址：http://srtm.csi.cgiar.org)。在该数据的官方网站上，将全球划分为若干个区域，每个区域为一个独立的文件，以“TIF”格式存储。文中使用的数据文件为“SRTM_59_08.tif”文件，该文件包含(109.9995838°—115.0004172°E,19.9995838°—25.0004171°N)范围内的所有区域，空间分辨率为0.00083°，在赤道地区约为90 m。我国珠江口附近的岛屿均包含在该文件中，文中特别截取了临近珠海东澳岛的区域进行试验。

使用的CFD软件为FLUENT通用计算流体力学软件，它的前处理器采用了GAMBIT，用于生成计算机辅助设计(CAD)模型并进一步划分网格，形成可被FLUENT读取的解域文件。具体技术路线如图1所示：1) 下载获取SRTM原始数据文件；2) 利用地理信息系统(GIS)工具读取原始文件，通过交互式方式选取并分割海岛地形高程数据，生成ASCII码格式的高程数据文本文件；3) 用FORTRAN语言编写GAMBIT脚本生成器，该生成器的功能主要是读取高程数据文本文件，并生成脚本，脚本可被GAMBIT读取自动执行；4) 利用GAMBIT读取生成的脚本文件并自动执行，形成高程点阵每个点阵为笛卡尔坐标系下坐标为(x, y, z)的点，然后通过交互式操作生成点阵的包络面，在建立包络面的基础上进一步生成体并划分网格；5) 输出.msh文件，该文件可被FLUENT等CFD软件识别和读取。

图1 SRTM数据处理技术路线

在生成网格文件之后，利用FLUENT进行风场模拟，其基本方程如式(1)—(3)所示：

(1)

(2)

(3)

2 结果与分析

2.1 地形数据与数据提取

利用GIS软件读取“SRTM_59_08.tif”文件，其可视化界面如图2a所示。由图可见，该文件包含了几乎整个广东省及南海的部分区域，北至南岭山地，南至雷州半岛，珠江口附近的小岛均被包含在文件中。通过逐级放大，可锁定目标小岛(如图2b中黑色方框所示)，裁取包含小岛的数据片段，即可获得生成CAD模型及相关脚本文件所必须的地形高程数据文件。

图2 锁定模拟区域及提取所需数据(a.地形文件；b.东澳岛)Fig. 2 Simulation area and the extracted data (a. Terrain file; b. Dong’ao island)

2.2 CAD模型与网格文件

利用FORTRAN程序编程得到GAMBIT脚本生成器。运行该生成器，可以根据上一步得到的数据文件生成一个脚本，脚本的每一行为一个自动操作命令，按照给定的坐标生成一个“点”。用于生成东澳岛CAD模型的脚本共有5 440行，在GAMBIT中运行脚本，相当于批处理执行生成“点”的操作，即建立了一个包含5 440个“点”的点阵，如图3a所示。进一步基于点阵，利用GAMBIT的面生成功能生成一个包络面，即可以较好地展现出海岛的地形特征(图3b)。

图3 GAMBIT生成的点阵(a)和基于点阵建立的包络面(b)Fig. 3 Dot matrix generated by GAMBIT (a) and envelope surface based on the matrix (b)

在包络面基础上，依次生成4个侧面和1个顶部的平面，并组合成一个“体”，即可作为模拟区域。顶部的高度一般应设置为地形高度的10倍，考虑到试验选取的东澳岛高度不到200 m，故解域的顶高设置为2 000 m。对比东澳岛CAD模型地面高程分布与卫星图片，可知文中建立得到的CAD模型无论从岛的海岸轮廓还是从地形起伏的特征来看，都较为一致。由此可见，尽管SRTM数据是开源数据，但基于该数据建立的海岛模型也很准确，足以为海岛风场的模拟研究奠定基础。

生成解域之后，须根据具体需求设置网格。文中使用了非均匀网格的设置，以图4a中的十字线为中心，分别向东西和南北两个方向逐渐拉伸网格距。东西和南北方向的网格数分别为200和160，水平方向上设置为30层，其中最底层的高度约为离地2 m左右。解域内下垫面的网格设置如图5所示。

图4 东澳岛CAD模型的地面高程分布(a)和东澳岛卫星遥感图像(b)Fig. 4 Ground elevation distribution of Dong'ao island by CAD model (a) and satellite remote sensing image of theisland (b)

图5 东澳岛风场模拟实验的地面网格设置

采用非均匀网格，是希望在东澳岛的东侧获得更高分辨率的数据，该侧为气流的来流方向，对于研究迎风侧的气流爬升等细节问题有帮助。而对于暂不关注的其他区域，采用较稀疏的网格设置，有利于节约计算资源，获得更快的模拟速度。

2.3 风场模拟结果

东澳岛东侧设有一个自动气象站，文中利用该站在2020年9月15日18:30左右监测到数据作为输入参数。该站在10 m高度处测得的2 min平均风速为5.7 m/s，入流风向角为69°(大致为东北偏东风向)。根据该观测数据，设入流风速的垂直廓线如式(4)所示：

(4)

式中，U(z)为z高度上的风速，U10为10 m高度的风速，z为距地面高度，0.2为指数风廓线在中性层结下的指数取值。将风速分别在东西和南北方向上分解，得到u分量和v分量，再将2个分量在垂直方向上的变化用FLUENT自带的用户自定义函数(UDF)编写出来。在正式模拟时，将事先编写的UDF载入模型，并在入流边界上赋值。

模拟分两步进行，第一步为稳态模拟，即不在模拟中加入时间偏导项。如前所述，在稳态模拟时使用雷诺平均的湍流闭合方案，realizable k-ε方案。待稳态模拟达到收敛以后，开启时间偏导项，并将湍流闭合方案设置为大涡模拟(LES)方案，积分3 600 s后，得到最终的模拟结果。

图6给出了模拟得到的距地面10 m高度的风矢量分布。由图6a可知，模拟结果较好的反映了复杂地形对风场的影响，在海拔较高的地方风速均较大，气流在爬越山峰的过程中，风速有明显增大的趋势，这与Li等(2014)之前在深圳排牙山地区的模拟结果类似。此外，图6a还表明，非均匀设置的网格有效地捕捉到了所关心区域的风场细节，在关注的东澳岛东侧地区，气流在接触岛屿时受地形阻碍损失动能形成的相对小风区较为明显。而在网格较为稀疏的其他地区，这一现象并不明显。图6b还给出了局部的模拟结果，由图可观察到更多的细节，特别是在背风侧，除了风速显著减小之外，风向也受地形影响出现了明显变化，在岛屿西侧的部分区域，风向已由东北偏东变为了偏北风，这是气流爬越和绕开岛屿后继续前进的表现。

图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进一步给出了模拟区域内沿主导风向垂直剖面上的风场分布，由图可见气流爬越岛屿的整个过程。在垂直方向上，风速随高度增高而增加的趋势十分明显。而在岛屿的山体背后，由于气流爬越山体损失部分动能，风速明显比山前要小。综上所述，利用文中提出的技术思路，可以较好地实现岛屿边界层风场的模拟，且模拟结果在物理上可以得到合理解释，这主要体现在风速的增强、衰减及风向的转变、扰动都符合预期，与地形的起伏有很好的对应关系，并与已有的研究结果相吻合。

图7 2020年9月15日18:30东澳岛风矢量(填色，单位：m/s)垂直分布模拟结果

4 结论与展望

文中提出了利用SRTM资料构建海岛的CFD计算模型的技术方法，并以珠江口的东澳岛为对象，进行了模型构建试验及边界层风场的数值模拟试验。实验结果表明：

1) SRTM资料作为开源的资料，可有效解析空间尺度在1 km内数量级的海岛地形；

2) 基于SRTM资料，可以构建用于CFD模拟的计算模型，并较为真实地刻画出岛屿的地形；

3) 利用计算模型进行的模拟试验表明，岛屿的边界层风场受岛屿地形影响十分明显，岛屿上空10 m高度的风场随地形增加而明显增加，在岛屿的背风侧，风力明显减弱，且风向有十分明显的变化。

必须指出，本研究仍然只是初步研究，重心放在了计算模型的构建上，对于岛屿边界层风场结构的科学分析比较粗浅。同时，由于目前观测资料的限制，模拟数据与观测数据的对比还未曾涉及，这些都是未来研究需要进一步弥补的地方。尽管如此，本文从技术上打通了从SRTM数据到CFD模拟的全流程，这些技术未来在海岛风资源评估、风电场选址、污染扩散等领域可以发挥作用。