沿海复杂地形台风登陆过程风场多尺度数值模拟

董浩天，陶 韬，杜晓庆

(1.上海大学，上海 200444；2.东京工业大学，日本 东京 152-8550)

0 引 言

我国东南沿海常年受到西北太平洋台风灾害的影响[1-2]。据国家应急管理部统计，仅2018年前三季度全国台风等风致灾害就造成4 712.2万人次受灾，203人死亡，直接经济损失854.2亿元。随着计算机硬件和软件的发展，近年来中尺度气象模式常用于台风过程的模拟和预测。WRF（Weather Research andForecasting model）开源中尺度气象模式采用了完全可压缩非静力平衡欧拉方程，并提供了大量物理参数化方案和数值方法选项，在大气边界层模拟和预报上明显优于传统大尺度气象模式[3]。Gu等[4]采用WRF模式模拟了2002年台风Rusa的中心路径和风速分布。Tse等[5]采用WRF模拟了台风过程近地风场，发现WRF可得到较好的平均风剖面和风速时程结果，但近地风速误差显着。Jiménez等[6]通过在控制方程中加入由地形复杂度定义的源项，提高了WRF模拟亚格子地形特征的能力。总的来说，WRF可以一定程度上预测/再现台风路径，并较好地模拟近地10 min平均风速时程；但由于控制方程形式、数值过滤、贴体网格和边界条件等原因，其对台风过程复杂地形近地湍流特征的模拟效果仍难以满足风工程的研究需求[7]。

复杂地形风场研究，特别是我国东南沿海地区复杂地形台风过程风场的研究，对灾害预防和风资源开发有着重要意义[8-10]。而同平坦地形相比，峡谷、高山、海峡等山区地形的特征尺度与大气边界层高度相近，易发生越山风、峡谷风和特征湍流等现象。随着湍流模型、网格技术、数值算法和计算机硬件的发展，计算流体动力学（CFD）方法已广泛应用于风工程领域并部分替代风洞试验，新版抗风规范[11]已对采用CFD虚拟风洞试验方法进行桥梁抗风研究给出了指南。Tamura[12]指出，大涡模拟方法（Large-eddy Simulation，LES）在非定常流动细节和计算资源消耗之间可达到较好的平衡，适用于结构抗风设计、复杂地形风场、城市风环境等风工程课题。Cao等[13]采用LES和浸入边界法模拟了湍流边界层风流经粗糙二维小山。复杂地形风场具备多尺度特征，即存在不同尺度、不同频率成分的大气流动。抗风规范[11]采用准定常假设，以10 min时距风速测量和规范风剖面、风谱构建平稳来流条件，可以较好描述小尺度流动现象，但忽视了桥址风场的多尺度特征。

为了克服中尺度模式在复杂地形近地湍流特征模拟方面的不足，Perivolaris等[14]学者提出了中尺度气象模式同小尺度CFD结合的多尺度数值模拟方法。如图1[15]所示，通过WRF嵌套网格求解大、中尺度三维风场低频信息；通过人工湍流生成技术补充高频风速成分，生成多尺度湍流风、温场，并作为边界条件输入CFD计算，求解近地风速时程。高频风速成分的再生成一直是多尺度数值模拟的难点。Wyszogrodzki等[16]直接将WRF结果经时间空间插值后输入LES城市地貌计算域，其上游地表起到了“数值粗糙元”的作用，湍流风剖面沿流向自然发展，但这一方法要求很长的上游计算域长度。Nakayama等[17]通过循环入流技术在WRF平均风剖面基础上进行人工湍流生成，降低了计算资源消耗。

图1 多尺度风场模拟技术路线Fig.1 A sketch of the multi-scale simulation of a wind field over a complex terrain

本文以2012年台风启德登陆广东东海岛为背景，采用多尺度风场模拟方法，研究了台风路径、近地风速时程和风剖面等，并同实测结果验证。研究内容可供东南沿海台风影响下的台风灾害预报、结构抗风设计和风力发电选址等参考。

1 台风模拟资料

1.1 台风启德概况

1213号台风启德（Kai-tak）于2012年8月13日在菲律宾以东洋面加强成为第13号热带风暴；8月15日21时加强为台风；8月17日4时30分前后在广东省湛江市登陆，登陆时中心附近最大风力达到13级；之后越过雷州半岛，进入北部湾并在越南再次登陆，最后于8月18日逐渐消散。图2显示了启德台风的中心路径。

图2 台风启德中心路径（来源：日本气象厅网站）Fig.2 The central route of typhoon Kai-tak 1213

1.2 气象观测资料

测站观测资料主要用于对近地风场模拟结果进行验证。广州省湛江市东海岛布置的测风塔记录了该次台风过境时的风速相关信息。东海岛紧邻雷州半岛，测风塔的位置为（21° 0' 44.90' N，110° 31' 33.88' E），该位置位于台风路径的北测。如图3所示，测风塔仪器安置高度分别在距离地表100 m、70 m、50 m、30 m高度处。测风设备为NRG #40风速计和NRG #200p风向传感器，输出数据包括风速和风向的10 min平均值、标准差、最大值和最小值。

图3 东海岛测风塔Fig.3 The wind observation tower at Donghai Island

1.3 气象再分析资料

气象再分析资料为多尺度数值模拟提供初值条件和边界条件更新。本次对台风启德模拟采用的再分析资料为美国国家环境预报中心（NECP）提供的FNL全球分析资料（Final Operational Global Analysis）。数据解析时间每6 h更新一次，格点分辨率为1°，包含了地表26个标准等力层（100～1 kPa）、地表边界层、对流层顶的要素信息。

1.4 地形地貌资料

在地表数据方面，对于WRF计算的外层d01和d02区域采用了美国地质勘探局（USGS）提供的1'高精度高程数据和地表覆盖数据；对于WRF内层的d03区域和LES计算区域，则采用国家基础地理信息中心提供的30 m精度ASTER GDEM地形数据和GLC地表覆盖数据。

2 多尺度模拟方法

2.1 WRF计算参数

中尺度数值计算采用Advanced Research WRF（ARW，即WRF用于科学研究的计算内核[3]）。ARW为基于全可压缩非大气静力平衡方程，采用贴地的压力定义的高度坐标系，模式顶的压力为5 kPa。WRF的计算域和网格参数如图4和表1所示。竖向共61层网格，其中在观测点位置最底部的10层高度分别为11.23、13.82、15.56、18.17、20.80、22.57、26.08、27.88、29.70、33.28 m。水平方向采用三层双向嵌套网格，即d01、d02和d03，最外层d01包括了台风启德登陆过程的区域，最内侧d03则主要覆盖了东海岛和周边地形。每层网格的格点数、间距和时间步长如表1所示，时间积分采用了Runge-Kutta方法。物理参数化方案方面，采用了WSM6微物理方案、Kain-Fritsch卷积云方案、Monin-Obukhov表面层方案、YSU大气边界层方案、RRTMG辐射方案和Unified Noah陆面层方案。

图4 WRF三层嵌套区域和台风启德中心路径结果Fig.4 Nested domains and typhoon routes

表1 启德模拟三层嵌套网格参数设置Table 1 Parameters in the WRF simulation

2.2 循环入流人工湍流生成技术

采用循环入流技术[18-19]进行高频风速成分的补充。如式（1）所示，空边界层数值风洞的入口速度分布(uin,vin,win)由WRF计算得到的平均速度剖面u0(z)和入口下游1倍边界层高度处循环断面速度(ure,vre,wre)的脉动部分组成，其中变量上划线代表时间平均值。

2.3 LES计算参数

小尺度LES风场计算的网格如图5所示。计算域的顶部采用对称边界，海拔高度统一为2 143 m。下部边界为无滑移壁面，海拔高度从0～78.28 m之间变化。北侧和南侧边界均采用对称边界条件，计算域南北方向宽800 m。东侧采用2.2节得到的速度时程作为速度入口条件，西侧则为出流条件，计算域东西方向长5 000 m。计算中主要模拟了来流台风边界层从东侧海平面向西越过沿海丘陵地区的过程，来流的平均方向为正西。水平方向网格划分，东西向Δx和南北向 Δz均为 Δx= Δz= 10 m，共 500 × 80个网格。竖向共划分150个网格，其中第一地表第一层网格高度Δymin= 0.143 m。六面体网格数量共计600万。

图5 小尺度LES计算网格示意图Fig.5 The computational domain and boundary conditions for LES

LES部分入口的平均速度采用WRF内层d03网格计算得到UTC时间2012年8月17日8：30附近的10 min平均速度剖面；风速脉动部分则通过2.3节的循环入流方法生成，入口风速剖面如图9和图10所示。计算中采用准定常假设，即入口平均风速分布不随时间变化。边界层顶的时均风速固定为Uref=27.82 m/s，边界层高度固定约为δ99= 600 m。在计算中通过改变空气黏性，以边界层高度计的雷诺数约为5.7 × 104。LES计算的时间步长固定为0.256 8 s，计算中克朗数保持在0.4以下。采样时间从t=1 027 s开始，此时以Uref= 27.82 m/s计的流体流经流场5.7遍；采样时长为600 s，对应流经遍数为3.3遍。

3 结果与讨论

3.1 WRF台风过程模拟结果

图4给出了WRF模拟得到的台风中心路径时程结果和日本气象厅提供的台风中心路径观测结果的对比。两者的时距均为6 h，从UTC时间2012年8月16日20：00至17日20：00，方向从东南到西北。受到初值偏差的影响[20]，d01初始段台风路径误差较大，达到60 km左右。经过网格双向嵌套以后，内层d03网格覆盖区域的台风路径精度较高，且明显优于外层d01和d02网格区域。

图6和图7给出了东海岛测站位置距离地表100 m高度处水平风速风向24 h时程的WRF结果与观测值对比。其中WRF和观测结果均为10 min平均值。时间按照UTC标准时间给出。WRF和观测得到的最大风速分别为30.3 m/s和33.6 m/s，数值较为接近，且均出现在17日12:30时附近。从17日11:00-12:00时范围的风速低谷区域以及对应的风向大致180度的变化来看，WRF可以较好地模拟台风眼过境影响下的风速、风向变化，但一定程度上高估了台风眼过境的时长。总的来说，WRF可以较好地模拟台风登陆过程的大气边界层风速时程，但模拟精度仍有待提高。

图6 东海岛测站100 m高度10 min平均风速Fig.6 10-min averaged wind speed at 100 m above the ground

图7 东海岛测站100 m高度10 min平均风向Fig.7 10-min averaged wind direction at 100 m above the ground

图8中蓝色点为LES计算得到的测站位置（x=0 m，z= 0 m）处地表10 m高度处x向风速和y向风速时程，而红色线为通过WRF平均风速+循环入流方法得到的入口（x= 2 500 m，z= 0 m）10 m高度处对应的风速时程。表2进一步给出了统计得到的10 min平均风速和湍流强度，其中湍流强度Iu=σu/Uref,Iv=σv/Uref，σu和σv分别为流向和竖向速度的标准差。可见，受到上游地形的影响，对于10 m高度处流向速度，测站处较入口平均速度下降8%，湍流强度则增加1倍左右；对于竖向速度，测站处出现0.73 m/s向下的平均速度，湍流强度稍弱于入口位置。图9给出了入口和测站10 m高度东西向风速的功率谱密度，可见测站位置的风速高频成分明显多于入口。

图8 测站位置和入口处10 m高度东西向风速u和南北向风速v时程Fig.8 Wind speed histories at 10 m above the ground

表2 测站位置和入口处10 m高度风速统计值Table 2 Wind speeds at 10 m above the ground at the inlet and observation tower

图9 测站位置和入口处10 m高度东西向风速功率谱密度Fig.9 The power spectra density of u at 10 m above the ground

3.2 LES复杂地形风场模拟结果

图10和图11给出了平均流向速度um和流向速度标准差σu的垂直剖面分布从上游到下游的变化。从地形来看，测站位置（x= 0 m）处于山顶偏背风侧。从入口（x= 2 500 m）开始，大约1 300 m左右的距离为海面。此外x= 1 000 m和x= -1 000 m两个位置则分别位于测站所处的山顶的上游和下游。从平均风剖面来看，测站处存在一定的越山加速现象，也即近地平均风速的提高；LES的结果比WRF得到的近地风速稍高，两者均高于观测值。从脉动风剖面来看，受到地表流动分离的影响，LES得到的近地脉动风速远高于WRF的结果，更接近观测值；但50 m以上的脉动风速下降明显。此外，观察到入口的脉动风速剖面在200～300 m段以及500～600 m段存在较大的值。导致LES模拟误差的主要原因可能有以下三点：1）WRF计算得到的平均风剖面误差传递到LES计算中；2）LES复杂地形风场计算和循环入流湍流生成段采用了相同且较低的雷诺数，与实际大气差异较大；3）公式（1）未对湍流向边界层顶的发展进行阻尼约束。

图10 上游（右）到下游（左）平均风速剖面变化（z = 0 m位置）Fig.10 Profiles of the mean wind speeds (z = 0 m)

图11 上游（右）到下游（左）脉动风速剖面变化（z = 0 m位置）Fig.11 Profiles of the fluctuating wind speeds (z = 0 m)

4 结 论

采用WRF中尺度气象模式和小尺度大涡模拟结合的复杂地形风场多尺度模拟方法，研究了2012年台风启德在广东东海岛登录过程的台风路径、100 m高度10 min平均风速风向时程、10 m高度近地风速时程、平均和脉动风剖面等结果，主要结论如下：

1）多尺度风场模拟方法可实现台风影响下沿海地区复杂地形风场时程的模拟，该方法不依赖规范风剖面假定或现场实测，对沿海地区台风灾害的预报和评估具有较大的现实意义。

2）台风登陆过程中受到沿海复杂地形的作用，近地顺流向湍流强度增加1倍左右，高频风速成分显著增加，同时出现了0.73 m/s向下的平均竖向风速。因此在台风登陆过程中，复杂地形局部可能产生远强于海面上的近地阵风，对抗风安全性产生不利影响。

3）后续研究可以从改善WRF风剖面模拟精度、提高雷诺数和修改循环湍流生成方法等三个方面入手，进一步提高近地风场的模拟精度。

致谢：感谢中国气象局、JMA、NECP、USGS、国家基础地理信息中心等机构提供的数据。