台风“山竹”影响下的城市风场数值模拟研究

2021-09-17廖孙策黄铭枫楼文娟肖志斌

空气动力学学报 2021年4期

廖孙策，黄铭枫，，*，楼文娟，林巍，肖志斌

(1.浙江大学建筑工程学院结构工程研究所，杭州 310058；2.浙江大学平衡建筑研究中心，杭州 310058)

0 引言

近年来我国台风灾害频发，根据中国气象局公布的数据，平均每年有6个以上台风级别的热带气旋影响我国，对我国沿海地区造成巨大的财产损失[1]。随着国民经济的不断发展，截至2020年末，我国常住人口城镇化率超过60%。对于沿海城市来说，台风带来的大风以及暴雨耦合独特的城市空气动力学行为，例如建筑群干扰效应[2]、热岛效应以及城市内部的复杂局部地形条件等，导致了严重的城市建筑群台风灾变破坏和各类次生灾害。台风作用下的城市风场特性也成为了国内外风工程研究的热点之一。

现场实测作为风工程研究的主要手段之一，被广泛的应用于台风风场特性的研究。不少学者[3-6]都开展了台风作用下城市风场特性的实测研究工作。近年来，数值模拟也日益成为风场研究的有效手段之一。Li等[7]利用Fluent软件与高精度城市建筑（Geographic information system，GIS）数据对晋江市某区域风场进行了CFD模拟，获得了详细、直观的风场信息，并与实测数据进行了对比。结果表明，结合GIS数据的模拟方法能够有效提高城市整体风场CFD模拟的效率；模拟结果与实测结果基本一致。但是，CFD模型受限于网格数量无法采用过大的计算域；此外，台风风场特性不同于良态风，仅依靠CFD模型无法准确获取台风作用下的入口边界条件。

不同于小尺度的CFD模型，WRF模式作为中尺度气象数值模式中应用最广泛的模型之一，具有开源、参数方案丰富等特点，能较好地预测、重现真实大气环境下的台风风场[8-10]。曾祥峰[11]利用WRF模式的三层固定嵌套网格对台风“约克”和“天鸽”进行了数值模拟，发现WRF模拟结果与实测结果相近，可以作为土木工程结构抗风设计的参考依据。此外，也有不同学者利用WRF耦合小尺度模式开展台风风场研究。黄铭枫等[12]基于WRF和CFD大涡模拟开展了杭州九堡大桥周边台风风场的数值模拟，并结合实测数据进行了验证。朱容宽与柯世堂[13]对厦门国际机场航站楼周边的台风风场进行了高时空分辨率WRF模拟，获得其边界层风速剖面作为航站楼风压CFD模拟的入口边界条件。为了细致考虑城市建筑的热力作用对WRF模拟结果的影响，CHEN等[14]提出了基于WRF的城市冠层模型（UCM），现有WRF版本已经集成该模块[15]。蒙伟光等[16]利用WRF耦合UCM模型，对广州地区的一次高温天气过程进行了数值模拟，结果表明该方法能准确模拟城市热岛的形成及分布，能改进对热岛强度的模拟结果。马敏劲等[17]利用WRF-UCM进行了降水过程模拟，发现UCM模型对水汽的分布和降水量都有一定程度的影响，耦合该模型后降水向城区靠拢。

已有研究表明，WRF模式能够有效重现台风的演化过程，捕捉台风风场特性。但是，上述关于WRF模拟台风风场的研究，还大多局限于对台风本身动力结构的分析或者是借助其他小尺度模式开展近地面风场模拟，此外对于WRF耦合UCM模型进行城市风场高分辨率模拟尚缺乏深入研究。考虑到台风风场结构的特殊性以及城市下垫面的复杂性，本文采用耦合UCM模型的WRF模式对台风“山竹”影响下的城市风场开展了多尺度数值模拟，并结合现场实测数据对城市台风风场模拟结果进行了验证。

1 台风“山竹”及现场实测概况

台风“山竹”（Mangkhut）于2018年9月7日20时在西北太平洋洋面上生成，并于9月11日升级为超强台风。2018年9月15日凌晨，台风“山竹”从菲律宾北部登陆后继续向我国南海移动并于9月16日17时在广东台山海晏镇登陆，登陆时中心附近最大风力14级，中心最低气压95.5 kPa。研究表明[3,18-19]，台风“山竹”是香港自1946年有风速记录以来最猛烈的台风，并存在明显的低空急流现象，给我国广东沿海地区造成了严重的经济损失与人员伤亡。

本文以深圳气象观测梯度塔（113.90°E，22.65°N）为中心模拟其周边区域在台风“山竹”过境时的城市风场。图1所示为深圳气象观测梯度塔，该塔由深圳市气象局和人居环境委共同建设，塔高356 m，是目前亚洲第一、世界第二的桅杆结构铁塔。图2展示了梯度塔周边地形及建筑分布情况。梯度塔位于深圳市宝安区铁岗水库，周边3～4 km范围内覆盖大片茂密的低矮树林，西侧及南北侧均有大片建筑群，东侧为海拔587 m的羊台山。梯度塔在10 m、20 m、40 m、50 m、80 m、100 m、150 m、160 m、200 m、250 m、300 m、320 m和350 m设置了共13个观测平台用于观测风速、气压、气温、降水等气象的梯度变化；其中10 m、40 m、160 m、320 m位置处设置有三维超声风速仪，采样频率为10 Hz，满足精度要求[20-21]。文献[20]与文献[21]利用该气象观测梯度塔分别开展了台风“天鸽”与台风“山竹”的实测研究。本文将结合文献[21]的观测数据对模拟结果进行验证分析。

图1 深圳气象观测梯度塔Fig.1 The meteorological gradient tower at Shenzhen

图2 梯度塔周边地形及建筑分布Fig.2 The terrain and buildings around the gradient tower

2 台风模拟试验设计

2.1 WRF模式及UCM模型介绍

WRF模式是由美国环境预报中心（NCEP）和多个科研机构共同开发的新一代中尺度气象数值预报系统。因其具有多重网格嵌套功能和丰富的大气物理过程参数化方案，该模式已经广泛应用于气象预报业务和相关学科的研究中[22]。

在WRF模式中参数化方案的选择对计算结果有重要影响[23,24]。WRF模式中集成了微物理方案、边界层方案、积云方案、辐射方案和陆面方案等多种参数化方案。对于城市风场模拟来说，大气湍流的求解尤为重要，针对网格尺寸较大的计算域，WRF模式采用行星边界层参数化方案来表征垂直方向湍流输送物理过程；而对于水平网格尺寸小于100 m的计算域，WRF模式支持大涡模拟（LES）模型，直接求解三维湍流中的大尺度涡部分，并通过亚格子应力模型替代小尺度涡，得到的结果更加接近真实大气环境中的湍流。

为精细化模拟台风“山竹”作用下的城市风场，本文在WRF模式的基础上耦合了UCM模型。不同于WRF中单一的城市土地利用类型，UCM模型将土地利用类型中的城市部分，根据建筑密度、道路宽度以及植被占比等细分为高密度住宅区、低密度住宅区和工业、商业用地三类[17]，并且用参数化的方式表征建造物的几何特征，同时也考虑了建筑物之间的遮挡、短波和长波辐射的反射、冠层的风廓线，以及屋顶、墙壁和路面的传热过程[15]，如图3所示。UCM模型对WRF中边界层的影响主要有热力作用与动力作用两种，均以通量的形式传递。热力作用方面，以感热通量为例，UCM模型分别计算植被与人为建造物的感热通量，再通过相应覆盖率加权计算得总感热通量。由UCM模型传递至WRF的总感热通量如下所示：

图3 城市冠层模型与竖向网格划分Fig.3 The urban-canopy modeling system and vertical mesh layer

式中QH为总感热通量；Fveg为城市中草地、树木等植被的覆盖率；Furb为城市中房屋、道路等人为建造物的覆盖率；QH,veg为植被的感热通量；QH,urb为人为建造物的感热通量，该项中的房屋屋顶、墙面和路面的感热通量通过给定参数分别计算[25]。总潜热通量与长波辐射通量的计算方法与总感热通量相同。在考虑城市冠层动力作用方面，UCM模型通过给定参数计算出房屋、道路、植被等的阻力系数和摩擦速度后，以动量形式传递至WRF边界层[15]。本次模拟中采用的UCM模型各项参数如表1所示，表中取值参考文献[15]与文献[17]。

表1 UCM模型参数表Table 1 Parameters in the UCM model

2.2 WRF计算设置

WRF模拟时间段为2018年9月13日8时至2018年9月17日8时，共120 h。为准确模拟台风“山竹”的动力演化过程，本文分两个步骤对台风“山竹”进行模拟。步骤一为台风演化过程模拟，如图4所示，步骤一共设置四层嵌套网格，前两层嵌套网格（Ad01、Ad02）为固定网格，包含了台风“山竹”在模拟时段内的完整移动路径；后两层嵌套网格（Ad03、Ad04）为移动网格，其会自动跟随台风中心移动以高精度求解台风近中心风场结构；四层嵌套域在垂直方向上均划分36层，第一层网格高度为26 m。步骤二为精细化城市风场模拟，如图5所示，以深圳气象观测梯度塔为中心建立5层嵌套固定网格（d03、d04、d05、d06、d07），最内一层嵌套网格（d07）水平分辨率达到49 m；垂直方向上采用非线性加密方案，第一层嵌套网格划分36层，后四层网格加密为72层，离地1 km内各层网格高度如表2所示。两个步骤之间通过WRF模式中的单向嵌套技术（ndown.exe）连接。步骤一中Ad02嵌套域的计算结果为步骤二的计算提供初始和每小时更新的气象边界条件。两个步骤嵌套网格的具体设置见表3。

表3 台风“山竹”城市风场模拟嵌套网格设置参数Table 3 The configuration of nested domains for the urban wind field simulation under the influence of typhoon Mangkhut

图4 台风演化过程模拟嵌套网格域Fig.4 Nested domains for the typhoon evolution simulation

图5 精细化城市风场模拟嵌套网格域Fig.5 Nested domains for the fine-scale urban wind field simulation

表2 1 km以下竖向网格分层高度Table 2 Vertical heights of each mesh layer within 1 km

合理的参数化方案选择对于精细化台风风场模拟至关重要。步骤一的台风演化过程模拟的微物理方案选择WSM6，行星边界层方案选择YSU方案，长短波辐射方案选择Dudhia和RRTM方案，陆面方案选择Noah Land Surface Model方案，对于计算域范围较大的两层嵌套域Ad01、Ad02选择Kain-Fritch积云方案。步骤二的精细化城市风场模拟的参数化方案中，微物理方案以及长短波辐射方案与步骤一保持一致，陆面方案选择Unified Noah LSM方案。在行星边界层方案的选择上，前三层嵌套网格（d03、d04、d05）选用YSU边界层方案，后两层嵌套网格（d06、d07）采用大涡模拟模型求解大气湍流。另外，在后三层嵌套网格（d05、d06、d07）中耦合UCM模型并设置相应对照组试验，研究UCM模型对城市风场模拟的影响。

除参数化方案外，WRF模式中的下垫面数据精度对模拟结果也存在一定影响。邹振操和邓院昌[26]在WRF模式中引入SRTM3和Asters高精度高程数据对香港地区进行了风场模拟，结果表明，与WRF原有数据相比，高精度数据在近地面风速模拟时能小幅提高模拟精度。本文在最内层的d05、d06、d07三层嵌套网格中采用Aster30 m精度地形高程数据和GLC2020-30 m精度土地利用类型数据，如图6与图7所示，相比于WRF原有数据，Aster和GLC数据能更准确地描述梯度塔周边的地形地貌及土地利用特征。

图6 d07计算域内地形高程数据对比Fig.6 The comparison of the topographical data in domain d07

图7 d07计算域内土地利用类型数据对比Fig.7 The comparison of land-usage categories in domain d07

此外，本文按照UCM模型的要求将土地利用类型中的城市细分为了高密度住宅区、低密度住宅区和工业、商业用地三类[17]，图8为广东沿海地区城市用地细分情况。

图8 广东及沿海地区城市用地分类Fig.8 The classification of urban land usage in Guangdong and coastal areas

3 模拟结果分析

3.1 台风演化过程模拟结果

图9所示为WRF模拟所得台风“山竹”移动路径与中国气象局（China Meteorological Administration，CMA）实测路径对比情况。模拟台风在前两天的移动路径几乎和实测路径一致，但在靠近菲律宾前移动方向产生误差，在擦过菲律宾后路径误差逐渐减小，其登陆位置几乎和历史登陆点一致。路径模拟结果表明，本次模拟采用的四层移动网格嵌套方案和选用的物理参数方案能够很好再现台风路径，大多数时刻的路径误差在50 km内。

图9 台风“山竹”模拟路径与CMA记录最佳路径数据对比Fig.9 The routes of typhoon Mangkhut from the simulation and CMA

图10给出了模拟所得台风“山竹”近中心最大风速、最低海平面气压与CMA实测结果对比情况。近中心风速对比结果表明，本次模拟准确地重现出了台风“山竹”在9月15日登陆菲律宾以及9月16日登陆我国广东省时的两次风速下降过程，模拟结果和实测结果相比总体风速大小保持一致，部分时刻风速略大于实测风速（图10（a））。由图10（b）可以看到模拟所得最低海平面气压和实测结果相比，上升幅度较小，但整体变化趋势一致，较好的模拟出了两次海平面气压的增大过程。综上所述，本文所采用的模拟方案准确地模拟出了台风“山竹”的移动路径和风场演化过程。

图10 台风“山竹”强度模拟结果与CMA数据对比Fig.10 Intensities of typhoon Mangkhut from the simulation and CMA

3.2 精细化城市风场模拟结果

通过五层嵌套计算域对深圳气象观测梯度塔周边风场开展了精细化模拟，最后两层嵌套域（d06、d07）在2018年9月16日17时（登录时刻）的10 m高度处风场如图11所示。d06和d07嵌套域使用了LES湍流模型并且耦合了UCM模型，可以看出，该方案精确模拟了深圳气象观测梯度塔周边的城市风场。如图11所示，台风“山竹”登陆时梯度塔周边最大风速已经超过30 m/s，由图可见，城市周边地形地貌以及建筑群的分布情况对局地风场有较大影响。如图11（a）所示，梯度塔周边山体导致了风速、风向空间分布不均匀，各个山顶位置处出现明显的高风速区，局部区域由于山体的遮挡出现低风速区，此外，部分位置的风向由于地形的干扰也出现了明显的转向现象。图11（b）为d07计算域得到的风场模拟结果，可以看到建筑稀疏区域的风速大于建筑密集区域；同时，从结果可以看到两块相邻建筑群之间形成狭道，气流从开阔地带穿过狭道时，过流断面缩小，根据流体连续性原理气流发生明显加速形成强风，这种现象被称为“狭道效应”；另外，由于来流方向上山体的遮挡，计算域下部位置出现了局地低风速现象。

图11 不同嵌套域10 m高度处水平风速云图（单位：m/s）Fig.11 Instantaneous contours of horizontal wind velocities (m/s) at 10 m above the groud in different domains

利用文献[21]中梯度气象塔分别位于10 m、40 m、160 m、320 m高度处的三维超声测风仪获取的10 min平均风速数据，与文献[17]中香港京士柏气象站（King’s Park，简称 KP）实测数据，对耦合 UCM 模型的WRF模式模拟城市风场的能力进行进一步的验证分析，同时比较UCM模型对WRF模拟结果的影响。数据时间段为2018年9月15日08时至9月17日08时，KP与梯度塔相对位置如图12所示。

图12 香港京士柏气象站地理位置Fig.12 The location of the King's Park weather station

图13（a-d）为WRF风速模拟结果与梯度塔实测结果对比情况，可以看到，WRF模拟结果与实测结果整体变化趋势一致，准确地模拟出了梯度塔位置处风速随着台风“山竹”的临近而不断增大，台风登陆后风速逐渐减小的过程。此外，从图中可以看出，WRF模拟所得风速稍大于实测结果，并且随着观测点高度的增加，模拟误差也有一定程度的增大，在160 m和320 m高度位置处模拟风速相较于实测风速减小较慢。产生误差的主要原因分析如下：台风中尺度模拟网格精度及地形解析精度的不足，会导致台风模拟风速偏大，这也是目前中尺度模式普遍存在的系统误差[27]；中尺度系统对于陆地表面复杂下垫面的解析精度不高，导致台风登陆后WRF模拟与实测之间的误差更加明显。此外，如图13（e）所示，由于梯度塔与KP相距约47 km，KP位置处最大风速出现时间与梯度塔相比稍有提前，梯度塔位置处的模拟风速也稍大于KP位置的风速。但风速、风向整体的变化趋势一致，这从一定程度上验证了模拟结果的可靠性。

图13 模拟风速与实测结果对比Fig.13 Comparisons of wind speeds between the simulation and observation

表4为不同高度位置处台风临近深圳气象梯度塔时（9月16日15时至16时）的最大风速比较情况，其中括号里的百分比数值表示相对误差。可以看到，单一WRF模式的计算结果与耦合了UCM模型均能较为准确地模拟出台风影响下梯度塔各个高度处的风速。在高度10 m和40 m的位置，耦合UCM模型后模拟风速更加接近实测风速，模拟精度有所提高；但随着高度的上升，二者之间的差距变小。这种现象的原因是UCM模型仅作用于前三层竖向网格（图3），所以高空位置处的风速模拟精度提升不大。另外，由于本次模拟采用的UCM模型较为简化，对深圳等复杂城市的真实特征描述不足。因此，计算结果与未耦合UCM的结果相比，风速模拟精度在低空略有提升，但总体相差不大。