高原气象

高原气象 >

2020 , Vol. 39 >Issue 1: 172 - 184

DOI: https://doi.org/10.7522/j.issn.1000-0534.2019.00005

论文

基于CFD降尺度的复杂地形风场数值模拟研究

  • 张嘉荣 , 1, 2 ,
  • 程雪玲 , 3
展开
  • 1. 宁波市海曙区气象灾害预警中心, 浙江 宁波 315153
  • 2. 南京信息工程大学大气物理学院, 江苏 南京 210044
  • 3. 中国科学院大气物理研究所, 北京 100029
程雪玲（1971 -）, 女, 辽宁大连人, 研究员, 主要从事大气边界层物理和大气湍流研究. E-mail:

张嘉荣（1993 -）, 男, 江苏泰州人, 硕士研究生, 主要从事风能资源数值模拟研究. E-mail:

收稿日期: 2018-07-07

  修回日期: 2019-01-10

  网络出版日期: 2020-02-28

基金资助

国家自然科学基金项目(41375018)

中国气象局全国风能资源详查和评价项目数值模拟专项

收起

Numerical Simulation of Wind Field in Complex Terrain based on CFD Downscaling

  • Jiarong ZHANG , 1, 2 ,
  • Xueling CHENG , 3
Expand
  • 1. Ningbo Haishu Meteorological Disaster Warning Center, Ningbo　 315153, Zhejiang, China
  • 2. School of Atmospheric Physics, Nanjing University of Information Science & Technology, Nanjing　210044, Jiangsu, China
  • 3. The Institute of Atmospheric Physics, Chinese Academy of Sciences, Beijing　 100029, China

Received date: 2018-07-07

  Revised date: 2019-01-10

  Online published: 2020-02-28

Fold

本文亮点

在计算机流体力学（Computational Fluid Dynamics, CFD）模式Fluent的基础之上修改传统的RANS湍流参数化方案, 选用大涡模拟（Large Eddy Simulation, LES）的湍流方案, 并引入温度层结的浮力模型, 得到修改后的Fluent模式。使用中尺度预报模式WRF和WRF耦合修改前后的Fluent模式模拟江西鄱阳湖地区内吉山测风站处2010年12月至2011年3月4个月的风速、 风向和2011年3月6日吉山站附近地形流场以及气温场的变化过程。结合观测资料的对比发现, WRF模式在耦合Fluent后, 对吉山单站点风速风向的模拟效果得到改善。WRF耦合Fluent的模拟方法在Fluent修改后（算例2）, 相比修改前（算例1）可以更精确地模拟吉山站各高度处的风速风向。4个月70 m高度平均风速WRF、 算例1和算例2的平均误差分别为3.963, 2.727和2.224 m·s-1。算例2与观测的70 m高度风向玫瑰图匹配程度总体上比算例1更高。算例2还可以模拟出白天不稳定状态下的热浮力效应和空间分布不均匀且较大的湍动能, 能模拟出夜间大气边界层的逆温现象, 而算例1则不能。以上结果均表明在对Fluent模式修改后, WRF耦合Fluent的模拟方法对吉山站附近风场的模拟效果得到改善。

本文引用格式

张嘉荣 , 程雪玲 . 基于CFD降尺度的复杂地形风场数值模拟研究[J]. 高原气象, 2020 , 39(1) : 172 -184 . DOI: 10.7522/j.issn.1000-0534.2019.00005

Highlights

The traditional RANS turbulence parameterization scheme is modified on the basis of the Computational Fluid Dynamics model (Fluent), and the large eddy model (LES) turbulence scheme is selected, then the buoyancy model of temperature stratification is introduced, finally the modified Fluent model is obtained.The Weather Research and Forecasting model (WRF) and WRF coupled Fluent before and after modified are used to simulate the wind speed and direction of the 4 months from December 2010 to March 2011 in Jishan wind station in Poyang Lake region, Jiangxi and the evolution of the flow field and the temperature field near Jishan station in March 6, 2011.According to the comparison of the observation data, it is found that: After coupling Fluent, the better WRF simulation result of the wind speed and wind direction in Jishan station is obtained.Compared with coupling original Fluent (Example 1), WRF coupling modified Fluent (Example 2) can provide more accurate simulation result of the wind speed and wind direction at each height in Jishan station.The average error of WRF, Example 1 and Example 2 for the 4 months 70 m height average wind speed are 3.963, 2.727 and 2.224 m·s-1 respectively.The Matching ratio to the observed wind direction rose of Example 2 is higher than Example 1 at 70 m height mostly.Example 2 not only can simulate the thermal buoyancy effect under the daytime unstable state and the spatial inhomogeneous and higher turbulent kinetic energy, but also can simulate the temperature inversion of the atmospheric boundary layer at night, while Example 1 can not.All the above results indicate that after the modification of Fluent model, the simulation result of the method that WRF coupled Fluent to the wind field near Jishan Station is improved.To sum up, this research provides a new method and basis for the prediction of the atmospheric flow field in the bottom of atmospheric boundary layer and the multiscale coupling simulation, and has also laid the foundation for the development and improvement of Fluent model in the future .

1 引言

真实的下垫面多具有复杂的地形起伏, 了解其对近地层流场结构的影响, 对于提高数值天气预报准确率、 研究大气污染问题、 气象灾害预警以及风资源详查等具有重要意义（李晓霞等, 2017; 曹杨等, 2017）。然而对于具有高度非均匀性复杂地形, 由于观测数据是一维的, 所能代表的范围非常有限, 若要获取二维或三维观测结果, 需要有大量实验设备（Taylor et al, 1987; Bechmann et al, 2009）, 因此类似的实验很少。除了现场实测之外, 获取复杂地形风场分布的方法目前主要还有风洞试验与数值模拟（肖仪清等, 2009）。风洞试验是通过建立复杂地形的缩尺模型, 在大气边界层风洞中测量风场的分布。综合考虑费用、 周期以及试验精度, 风洞试验较现场实测方法并不具明显优势, 只适用于典型地形的风场验证研究（肖仪清等, 2009）。而数值模拟方法并不受上述条件限制, 利用这一手段获取复杂地形条件下高分辨率的近地层流场资料, 显得尤为重要。
国内外针对复杂地形条件下近地层流场结构的研究可包括诊断模式和预报模式两类。诊断模式与预报模式两者在方程是否使用时间项上存在差异。诊断模式的方程中不含有时间项, 它以实测数据出发, 并通过计算得到具有稳定状态的流场。预报模式的方程中含有时间项, 它以初始场开始, 考虑边界条件的影响, 通过计算得到气象场的时间演变。
诊断模式可以分为两种, 包括线性模式和质量守恒风场调整模式。线性模式通过线性化方法简化Navier-Stokes方程后计算得到定常解, 而质量守恒风场调整模式是通过质量守恒定律计算分析各种气象场数据。目前比较常见的诊断模式包括WAsP(Wind Atlas Analysis and Application Program)（Mortensen et al, 1993）、 MSFD（Mixed Spectral Finite Difference Model）（Beljaars et al, 1987）、 CALMET（California Meteorological Model）（Wang et al, 2008）等。
预报模式则是通过各类气象因子的守恒方程在时空尺度上的积分来预报气象场的时空演变过程, 守恒方程包括质量方程、 动量方程、 热量方程和标量方程等。不同的预报模式采用的近似方案和对物理过程的参数化方案有所不同, 各种近似方案如不可压缩、 静力平衡、 Boussinesq假设等, 参数化方案如湍流、 通量、 边界层等。它们都有各自的局限性（Thunis et al, 1996）。目前被广泛使用的预报模式有WRF（Weather Research and Forecasting）、 MM5（Mesoscale Model version 5）、 RAMS（Region Atmosphere Model System）、 ARPS（Advanced Regional Prediction）等, 其边界层、 近地层和次网格参数化方案对近地面大气流场的模拟结果影响较大。
仅使用预报模式, 很难满足复杂地形条件下近地层流场的模拟精度要求。即使模拟精度可以达标, 对于区域性模拟需耗费大量机时。所以, 需要结合预报模式与诊断模式两类模式, 即使用动力学降尺度的手段来实现近地层风场高精度模拟。例如美国国家可再生能源实验室（U.S, National Renewable Energy Laboratory, 2006）用中尺度模式MASS与小尺度模式Mes-Map结合的模式系统完成了美国、 巴西、 中国东部等国家和地区水平分辨率1 km×1 km的风能资源模拟; 加拿大气象局将中尺度模式MC2（Mesoscale Compressible Community model）与小尺度模式Ms-micro相结合建立了精细化风能资源数值模拟评估模型WEST(Wind Energy Simulation Toolkit), 并绘制了5 km×5 km分辨率的加拿大全国风能资源图谱（Yu et al, 2006）, 何晓凤等（2010）采用MM5模式与WindSim结合进行了都阳湖区域风能资源的数值模拟研究等, 周荣卫等（2018）采用中尺度气象模式WRF结合牛顿松弛逼近Nudging资料同化技术, 实现哈密地区水平分辨率1 km的近地层风场数值模拟计算。刘郁珏等（2018）利用WRF模式与LES结合实现了对2022年北京冬奥会小海坨山赛区边界层风场的精细数值模拟。
近年来, 随着计算机能力的增强, 以往用于空气动力学精细流场计算的CFD（Computational Fluid Dynamics）模式越来越多的在气象领域得到应用, 很多学者将其应用于复杂地形近地层风场的模拟中（何晓凤等, 2010; 程雪玲等, 2005, 2006 Stangroom, 2004; Benchmann, 2006; 李磊等, 2010; 方艳莹等, 2012; 刘丽珺等, 2018; 姜平等, 2019）。在此基础上, 人们开始研究用中尺度预报模式和CFD模式结合进行复杂地形风场的数值模拟（何晓凤等, 2010; 方艳莹等, 2012）。
Fluent模式是国际知名的CFD模式, 它采用以有限体积法(Finite Volume Method)为核心的数值计算方法, 具有较高的收敛精度, 可用于结构化和非结构化网格方案的计算求解。Fluent模式的前处理模块Gambit有着很强大的网格划分和建模功能, 可以生成模拟区域内各种复杂的贴体网格。同时, Fluent还集成了非常丰富的参数化方案以及各种物理模型, 以应对各种不同流体流动问题的模拟。另外, Fluent还提供了用户自定义函数(User Defined Functions, UDF)接口, 用于链接外部数据和增加物理化学模型。因此, 本文将使用Fluent模式来模拟复杂地形条件下的近地层风场。
Fluent的湍流模式既有多种基于RANS方程的湍流闭合方案, 也有基于LES模式的各种亚格子(sub-grid)模型。Wilcox（1998）认为, RANS方程的湍流闭合方案是基于Boussinesq假设, 假定雷诺应力只与当地的平均风速有关, 且湍流是各向同性的。而这一假设在大气边界层中是不准确的, 尽管二阶k-ε湍流闭合方案能够描述大气湍流的各向异性行为, 也能计算应力输运项, 但耗散方程（ε方程）中仍需人为设置参数, 是造成计算误差的一个很大的来源。Kim et al（2000）采用基于RANS方程的不同的湍流闭合方案模拟了Askervein山的风场, 认为RNG k-ε湍流闭合方案能够精确地计算出风场的平均速度、 气流的分离和再入, 以及湍流特征量。Summer et al（2010）认为基于RANS方程的数值模拟方法虽然能够准确地计算风场的平均值, 但对于湍流的脉动值却并不准确, 尽管采用RNG k-ε湍流闭合方案能描述气流的分离再入, 却不能得到湍流的间歇结构、 涡旋的产生和输送过程, 这是RANS方程无法描述脉动运动导致的。因此, 更倾向采用LES模式进行复杂地形的流场计算。LES模式经亚格子过滤后的小尺度湍流更接近各向同性, 此时采用湍流闭合方案计算更合理和准确（Wood, 2000; Chow et al, 2009）。
CFD模式被用于大气边界层复杂地形模拟之后, 在发挥其精细流场模拟优势的同时, 人们越来越关注其对边界层温度层结的模拟。由于大气的温度、 密度、 水汽等性质在垂直方向有明显变化, 所以大气是一种层结流体, 空气的垂直运动和大气层结的不同类型密切相关。而空气运动受复杂地形的影响, 使得在平坦下垫面成立的相似理论等不再能够使用（Coppin et al, 1994; Hunt et al, 1991; Weng, 1997）。随着CFD模式的日渐成熟, 一些CFD模式也增加了对温度层结的模拟（Meissner et al, 2011; Koblitz et al, 2011）。Fluent采用Boussinesq密度假设求解非定常可压流动, 可以通过用户自定义函数(UDF)接口, 增加温度层结的物理模型。如Kristof et al（2007）通过增加动量方程、 能量方程和湍流输运方程中的源项, 模拟了城市热岛现象, Russell（2009）和Leblebici et al（2011）用Fluent模拟了复杂地形不同稳定层结的边界层流场等。
江西省的鄱阳湖地区以山地丘陵为主, 本文以江西省鄱阳湖地区吉山塔附近区域为模拟对象, 并根据模拟方案中Fluent模式湍流方案选择的不同和温度层结浮力模型的引入与否, 分别对模拟对象实施模拟, 旨在改进现有的Fluent模式, 探索出一种能够更准确模拟复杂地形条件下近地层流场的WRF耦合Fluent模式模拟方法。

2 资料选取与方法介绍

2.1　模拟区域与观测资料

模拟区域如图1[该图基于江西省自然资源厅标准地图服务网站下载的审图号为赣S(2017)043号的标准地图制作, 底图无修改]所示, 鄱阳湖(115°49′E -116°46′E, 28°24′N -29°46′N)位于我国江西省北部, 长江中下游南岸。鄱阳湖地区属亚热带湿润季风气候, 降水充足, 气候温和。鄱阳湖平原由赣、 抚、 信、 饶、 修等河流冲积而成, 地势平缓, 海拔多在50 m以下, 边缘部分有相对高度20~30 m的红土坡状地。全湖现有岛屿共41个。
图1 模拟区域和测风塔的位置

Fig.1 The simulation area and the location of the windmeasurement tower

在这片区域设有多座风能测风塔, 永修县吉山塔(116.07152°E, 29.23060°N)位于鄱阳湖西岸的吉山岛上永修县, 海拔为89 m, 该区域为浅山丘陵地貌, 面积约17.5 km2。测风塔高度为70 m, 观测层数共计4层。目前获得吉山站的资料为2010年12月至2011年3月期间的观测资料（风速和风向）。

2.2　模式介绍

2.2.1　 WRF模式

WRF模式是新一代高分辨的中尺度天气预报研究系统, 是由美国国家大气研究中心(NCAR)、 美国国家环境预测中心（NCEP）等多个研究部门及美国大学的科学家共同开发建立的。其中WRF-ARW（Advanced Research WRF）适合用来研究尺度跨度数百米到数千公里的问题, 包括数据同化研究、 预报研究、 区域气候研究等等。本文采用的是WRF3.7版, 它的介绍可由以下网址获得（http: //www2.mmm.ucar.edu/wrf/users/docs/user_guide_V3.7/contents.html）。

2.2.2　 Fluent

2.2.2.1　基本控制方程

Fluent使用的控制方程都是由Navier-Stokes方程演变而来。在模拟小尺度的流场时, Fluent假设空气运动是3D的不可压缩稳定流动, 忽略地转偏向力（科氏力）的影响。
将空气运动速度ui和气压P分解为平均值和脉动值, 平均值用上标“—”表示, 脉动值用上标“′”表示, 即:
U i = u i ¯ + u i ',
P i = p i ¯ + p i ',
Navier-Stokes方程可以重新写为:
∂ ∂ x i ( ρ u i ¯ ) = 0,
∂ ( ρ u i u j ¯ ) ∂ x j = - ∂ p ∂ x i + ∂ ∂ x i μ ∂ u i ¯ ∂ x j + ∂ u j ¯ ∂ x i - 2 3 δ i j ∂ u l ¯ ∂ x l                 + ∂ ( - ρ u i ' u j ' ¯ ) ∂ x j + ρ g ¯  ,
式（3）和（4）就是雷诺平均Navier-Stokes方程（RANS）, 其中, ρ为空气密度; δ i j为克罗内克（Kronecker）张量; g为重力加速度。式（4）中多了一个未知项雷诺应力张量 - ρ u i ' u j ' ¯, 它代表着速度脉动量的相关性, 该项使RANS方程不闭合, 需要将其参数化, 使方程闭合。
Fluent集成了多种湍流雷诺应力参数化方案, 包括k-ε方案, k-ω方案, 雷诺应力输送模型以及大涡模型（LES）。其中k-ε方案对计算精度和计算资源进行合适的匹配, 确保可靠的结果和稳定, 是最为成熟的做法。而LES方案能够提供最好的参数化结果, 不过比较耗计算资源。
在三种k-ε方案[标准k-ε方案（Launder et al, 1972）、 RNG k-ε方案（Yakhot et al, 1986）和realizable k-ε方案（Shih et al, 1995）]中, Li et al（2006）发现realizable k-ε方案模拟复杂地形的大气运动的效果最好。
大涡模型主要有两种亚格子模式（Smagorinsky模式和Lagrange动力模式）。在近地层大气流场的LES数值模拟中证实Lagrange动力模式优于Smagorinsky模式, 而且计算量增加不多（约15%）（Tseng et al, 2006; Shi et al, 2008）, 因此选用Lagrange动力模式。在两种lagrange动力模式（WALE模型和Kinetic-Energy Transport模型)中, 选用Kinetic-Energy Transport模型。
在后面的算例中将分别使用realizable k-ε方案和LES模型两种湍流方案进行模拟, 并对这两种方案进行比较。

2.2.2.2　Fluent静态数据

Fluent所需静态数据为建网格所采用的高精度地形高程数据, 采用美国航天局（NASA）与日本经济产业省（METI）于2003年共同推出了最新的先进星载热发射和反射辐射仪全球数字高程模型ASTER（The Advanced Spaceborne Thermal Emission and Reflection Radiometer）GDEM（Global Digital Elevation Model ）。该产品每1度经纬度方格划分一个文件, 30 m（约1弧秒）精度, 投影系统为WGS84（World Geodetic System-1984 Coordinate System）的GeoTiff地球经纬度输出格式。可以从NASA网站（https: //wist.echo.nasa.gov/api/）公开下载, 也可以登录METI的地面数据系统查询下载（http: //www.meti.go.jp/）。

2.2.3　 WRF驱动Fluent流程

WRF驱动Fluent的主要思想是将WRF相关网格点的输出气象要素（主要是风速、 风向以及气温）通过插值当成Fluent网格的边界条件, 驱动Fluent求解局地风场。在驱动过程中, Fluent每隔1 h或0.5 h更新一次边界条件, 并且假设这段时间内Fluent的外场不变, 初始场为上一时刻的流场。这是一种单向的驱动。

2.3　模拟设置

针对如图1所示鄱阳湖核心区域和吉山塔的位置分别建立WRF和Fluent的网格, 并设置好相关参数。

2.3.1　 WRF设置

WRF对模拟区域采用了三层网格嵌套。最外层区域格距为54 km, 格点数为68×68, 第二层区域格距为18 km, 格点数为91×91, 范围覆盖了东部南部大部分地区, 第三层区域格距为6 km, 格点数为82×64, 范围主要覆鄱阳湖地区。
模式使用的地形资料为30”×30”的USGS格点资料, 其分辨率不超过0.9 km。模式在垂直方向上使用sigma坐标, 共计32层, 按照上疏下密的方式划分, 其中1 km高度往下共计12层。模式最顶端气压为50 hPa（高度约为21 km）。初始场为NCEP再分析格点资料, 设置积分步长为20 s。物理方案的选择依次为Lin skeme（微物理过程方案）, Monin-Obukhov (Janjic Eta) Scheme（近地层方案）, Unified Noah land-surface model（陆面过程方案）, MYJ scheme（边界层方案）, Rrtm skeme（长波方案）, Goddard short scheme（短波方案）, NONE（城市冠层方案）, Kain-Fritsch scheme [积云方案（第三层不采用积云方案）]。

2.3.2　 Fluent设置

根据图2, 截取吉山塔周围3.6 km×3.6 km的GDEM, 构建7.5~200 m水平分辨率的网格, 设置边界条件, 将cas文件读入Fluent后设置好各项参数。
图2 Fluent接受WRF边界条件求解更细微尺度风场的技术路线

Fig.2 The technical route that Fluent accepts the WRF boundary condition to solve the more subtle wind field

2.3.2.1　准备地形数据

要生成Fluent模拟计算近地层大气流场所需要的高分辨率网格, 需要把模拟区域的30 m分辨率的ASTER GDEM数据准备好。首先, 要对数据的准确性进行验证, 将已知经纬度地点的海拔和数据集内对应点的高度作对比, 或者通过对照数据集内的显著地形（如山峰或洼地）来验证。然后, 需要对数据进行正交转换。因为从NASA或METI下载的ASTER GDEM数据是WGS84地球基准的地理坐标系统, 为了方便后续建立网格, 将经纬度坐标投影转化成通用横轴麦卡托投影坐标（UTM, Universal Transverse Mercator）(单位: m), 不过由于投影方式的影响, 转为的UTM坐标地形数据不再为正交坐标系, 需通过matlab处理使之正交。之后, 使用matlab截取模拟区域范围内的地形数据。
图3(a)表示的是截取的鄱阳湖吉山附近地形, 图中黑线为吉山塔的位置, 以吉山塔为中心四周范围为3.6 km×3.6 km。中心区域为一座坡度较小的小山（海拔100 m左右）, 四周为平坦地形（海拔5 m左右）。
图3 模拟区域

Fig.3 The simulation area

2.3.2.2　产生网格

按照程雪玲等（2006）所述方法, 使用perl软件编程将截取的地形数据写为相应的Journal脚本, 并且输入到前处理器Gambit中, 脚本执行建点与建面的命令, 将三个控制点连接形成一个控制面, 然后将所有的控制面相互连接在一起即形成了地形虚面。通常为了避免侧向边界压力场出现异常, 在模拟区域的四周再延展一定距离的平坦面[图3(b)], 延展后的网格大小为6 km×6 km。
计算网格可分为结构性和非结构性网格两类, 结构性网格便捷易于建立, 计算效率高, 不易出错。本文生成的地形曲面是正方形, 因此生成结构性网格用于计算。首先在地形虚面上建立面网格, 设置水平方向分辨率为50 m×50 m, 其中算例2在中心1.2 km×1.2 km的区域使用分辨率为7.5 m的加密网格, 以此充分描述湍流的过程。垂直方向上的网格划分算例1、 2相同, 第一层网络高度为1 m, 1 m外网格逐渐拉升, 最顶层网格高度为1 km, 最后生成六面体网格, 如图3(b)所示（仅显示算例2的网格）。总网格数算例1为84×84×50, 算例2为220×220×50。

2.3.2.3　边界条件

假设空气运动是3D的不可压缩稳定流动, 忽略地转偏向力（科氏力）的影响。设定空气密度为常数, 为1.18 kg·m-3), 动力粘性系数为1.7894×10-5 N·s·m-2), 地表粗糙度为0.01 m。
针对根据实际地形建立的六面体结构性网格有东西南北上下共计6个面需要设定边界条件。通过Fluent的UDF接口来提供入口边界条件。
根据WRF最内层网格输出, 插值确定吉山模拟区域侧向边界格点的风速、 风向及气温, 每1 h更新一次。其中算例1使用的边界条件为风向风速, 算例2使用的边界条件为风向风速和气温, 算例2模拟区域的地表温度由地表热平衡方程计算:
∂ T s ∂ t = C T ( R n - H - L E ) - 2 π τ ( T s - T 2 ),
式中: T s为地面温度; T 2为深层土壤温度; τ是根据土壤条件确定的松弛时间; C T是地表的热系数。
根据风向设置来流的方向为速度入口（velocity inlet ）边界, 流出方向为出口（outflow）边界。上边界为对称（symmetry）边界, 下边界无滑移壁面（no-slip wall）。

2.3.2.4　算例设计

设计的算例包括WRF直接模拟和WRF耦合Fluent模式的两个算例（算例1、 2）共计三种算例, 其中算例1、 2在WRF模拟后耦合了Fluent, 在Fluent湍流模型的选择上, 算例1选择realizable k-ε模型, 算例2选择大涡模型（LES）。
Fluent中LES湍流模型的动量方程本身没有浮力项[其动量方程如式（6）所示（详情可参考Fluent 6.3 User’s Guide）], 无法分析浮力效应对大气流动的影响, 需在方程右侧人为加入浮力项 θ ¯ - θ 0 ¯ θ 0 ¯ g δ i 3[如式（7）所示]。本文先根据式（7）编写好该动量方程的UDF（用户自定义函数）的C语言程序, 然后通过Fluent的UDF接口链接, 以此实现对浮力模型的引入。此外, 算例2中的气温通过LES模型的温度对流扩散方程[式（8）]和亚格子模式Kinetic-Energy Transport模型来求解。
∂ u i ¯ ∂ t + u j ¯ ∂ u i ¯ ∂ x j = - 1 ρ ∂ p ¯ ∂ x i + ν ∂ 2 u i ¯ ∂ x j ∂ x j + ∂ τ i j ∂ x j + f i ¯,
∂ u i ¯ ∂ t + u j ¯ ∂ u i ¯ ∂ x j = - 1 ρ ∂ p ¯ ∂ x i + ν ∂ 2 u i ¯ ∂ x j ∂ x j + ∂ τ i j ∂ x j +   θ ¯ - θ 0 ¯ θ 0 ¯ g δ i 3 + f i ¯    ,                                                  （ 7 ）
∂ θ ¯ ∂ t + u j ¯ ∂ θ ¯ ∂ x j = κ ∂ 2 θ ¯ ∂ x j ∂ x j + ∂ τ θ j ∂ x j + S θ ¯,
式中: p ¯和 θ ¯分别为可解尺度的气压和温度; κ为热扩散系数; S θ ¯为热源/汇项; f i ¯为除浮力外的质量力; θ ¯ - θ 0 ¯ θ 0 ¯ g δ i 3为浮力项; θ 0 ¯为参考温度; τ i j = u i ¯ u j ¯ - u i u j ¯为亚网格应力; τ θ j = θ ¯ u j ¯ - θ u j ¯为亚网格热通量。本文设计的三个算例的详细情况见表1。
表1 三个算例详细情况

Table 1 The details of three examples

算例 WRF 算例1 算例2
是否耦合Fluent 否 是 是
是否引入浮力模型 - 否 是
Fluent中湍流模型 - realizable k-ε model LES model

3 结果分析

3.1　一次大风过程的模拟

下面模拟2010年12月13日08:00(北京时, 下同)至16日08:00时段的天气过程, 定量地分析WRF模式直接模拟和算例1、 2模拟鄱阳湖吉山站的风速与观测资料的对比。
根据2010年12月13 -16日这4天的天气记录, 发现这段时间内出现了一次全国范围的降温寒潮过程。从图4可以看出, 鄱阳湖地区的这次大风过程, 从0 h（13日08:00）一直到54 h（15日14:00）间风速不断地增大, 于15日14:00风速达到最大值为23 m·s-1, 随后风速迅速减小。
图4 2010年12月13日08:00至16日08:00吉山70 m高度处风速模拟结果与观测资料对比

Fig.4 The comparison of wind speed simulation data andobservation data at 70 m height in Jishan station from08:00 on 13 to 08:00 on 16 December 2010

结合图4和表2可以看出, 相对WRF模式直接模拟, 算例1能够较大改善平均误差, 相关系数也略有改善, 而均方根误差则稍变大。算例2在三个指标上都优于WRF模式直接模拟和算例1。这说明, 算例2对大风过程的模拟优于WRF模式直接模拟和算例1。
表2 2010年12月13日08:00至16日08:00吉山70 m高度处风速模拟结果的各项评估指标

Table 2 The evaluation indexes for the wind speed simulationresults from 08:00 on 13 to 08:00 on 16 December2010 at 70 m height at Jishan station

WRF 算例1 算例2
平均误差BIAS/(m·s-1) 4.602 3.755 3.104
均方根误差RAMS/(m·s-1) 3.956 4.133 3.423
相关系数R 0.785 0.791 0.883

3.2　鄱阳湖模拟总体结果与分析

总体模拟时间段为2010年12月至2011年3月, 其中WRF模式每7天（共168 h）模拟一次, 其中每次第一天的模拟结果不采用。将WRF、 算例1和算例2的模拟结果分别与观测资料进行对比。模式结果和观测资料都是1 h平均后的结果。

3.2.1　总体结果

图5为吉山站70 m高度处的模拟风速与观测风速的时间序列图, 为了显示清楚, 这两张图显示了1440 h的对比结果。表3为对三种不同的模拟结果进行进一步的统计分析。
图5 吉山70 m高度处风速模拟结果与观测资料对比

Fig.5 The comparison of wind speed simulation and observation data at 70 m height at Jishan station

表3 吉山70 m高度处风速总体模拟结果的各项评估指标

Table 3 The evaluation indexes for overall simulation ofwind speed at 70 m height at Jishan station

WRF 算例1 算例2
平均误差BIAS/(m·s-1) 3.963 2.727 2.224
均方根误差RAMS/(m·s-1) 3.708 3.351 2.752
相关系数R 0.710 0.682 0.778
结合图5与表3比较可知, WRF模式直接模拟的风速偏高, 相比之下, 算例1模拟风速平均误差有显著改善, 与此同时, 均方根误差也得到改善, 相关系数略有下降。而算例2在平均误差、 均方根误差和相关系数不仅相对算例1都有改善, 相关系数也大于WRF模式直接模拟。由此可知算例2的模拟效果比WRF模式直接模拟和算例1要好。

3.2.2　月统计结果

3.2.2.1　月平均风速

分别计算2010年12月至2011年3月内每月的吉山站6层高度（4, 10, 30, 50, 70和100 m）上的平均观测风速、 WRF模式输出（吉山对应格点）平均风速、 算例1输出（吉山对应格点）平均风速和算例2输出（吉山对应格点）平均风速, 给出相对误差。WRF模式直接输出和算例1、 2输出（吉山对应格点）的相对误差见图6。
图6 吉山站6层高度月平均风速三个算例模拟的相对误差

Fig.6 The relative error for monthly mean wind speed simulation result of three example at six height at Jishan station

从图6可以看到, WRF模拟的月平均风速在各层高度都偏高, 其中4 m和10 m模拟的相对误差较大。经过CFD模式处理之后, 吉山站6层高度月平均风速模拟的相对误差明显减小, 尤其是算例2, 相对误差得到显著改善。6层高度4个月平均风速算例1的相对误差分别为30.9%, 24.7%, 14.7%, 16.5%, 16.3%和10.2%; 算例2的相对误差分别为26.9%, 22.5%, 12.8%, 12.3%, 13.3%和9.3%。
图6显示的是数值模式直接输出的结果（格点的最邻近插值）, 未经过后续的订正。

3.2.2.2　风向玫瑰图（70 m）

由吉山站70 m高度处风向玫瑰图（图7）可知, 2010年12月至2011年3月吉山站风向以北东北方向为主, 2011年1月风向频次最多, 达到了70%, 算例1和算例2对风向的模拟较好, 不过二者对2010年12月份非盛行风风向的模拟并不理想。经对比可以看出, 算例2输出的2011年1月和3月两个月风向频次比要算例1更接近观测数据, 其余两个月两者基本相似。
图7 吉山站70 m高度处观测和算例1、 2输出的风向玫瑰图

Fig.7 The wind direction rose diagram at the 70 m height at Jishan station: observation, example 1 and 2

3.3　典型时刻的流场和大气温度场时空特征分析

在与观测数据对比评估的基础之上, 通过进一步研究吉山站附近地区内流场和大气温度场的时空特征, 评估算例1和算例2两种模拟方式的模拟效果。选择2011年3月6日14:00和24:00两个典型时刻, 分析1 h平均的流场和温度场的演变特征。
图8(a)、 (b)分别显示14:00计算域侧向边界处WRF模拟的风速风向和温度剖面, 图8(c)、 (d)是域内算例2模拟此时的风速和温度场, 其中XZ界面位于Y=1800 m处。从图8(a)中可看出, 此时域内的风向为东南风, 图8(c)显示在距地面600 m处有较高的风速, 在山后（按主流风向区分山前山后）尤为显著, 可达15 m·s-1。山后近地面部分区域气流风速显著减小, 近地层上空有明显的垂直对流, 形成具有一定高度的垂直环流。从图8(d)中可以看出, 地面此时受强烈太阳辐射的作用, 温度升高, 而气温随着高度的增加而降低, 在300 m左右以上的高度, 温度的变化很小, 由此可见日间不稳定近地层空气在垂直方向上的混合作用。近地面空气在山后增温, 这是因为气流过山后风速减慢, 受到地表感热的影响, 温度升高。山后近地面上空温度隆起的区域与风速中出现的垂直环流表明垂直对流在浮力的作用下得到充分发展, 这体现了浮力效应对大气流动的影响。
图8 14:00风速风向(a)、 气温(b)的侧边界条件与算例2(c, d)和算例1(e)的模拟结果

Fig.8 The lateral boundary conditions of wind (a) and air temperature (b) and the simulation resultsof example 2 (c)， (d) and 1 (e) at 14:00

图8(e)显示的是域内算例1模拟14:00的风速（此时Fluent不考虑热对流对大气流动的影响, 所以不能提供温度场的分布）, 从图8(e)中可以看出, 在垂直方向上, 气流遇山后受迫抬升, 风速在山后有所减小。山后近地面上空气流较平直, 并没有出现明显的垂直环流。说明该算例不能模拟出白天气流过山后减速增温形成的热对流。
图9(a)、 (b)分别显示24:00计算域侧向边界处WRF模拟的风速风向和温度剖面, 图9(c)、 (d)是域内算例2模拟此时的风速和温度场。从图9(a)中可以看出, 此时的主流风向为西南风, 与14:00的东南风发生明显的改变, 此时的风速比14:00略有减小, 山后没有出现14:00那样的垂直环流, 山后近地层上空气流较为平直[图9(c)]。从图9(d)中可以看出, 由于夜间太阳辐射消失, 地面长波辐射使得地表降温, 在700 m高度以下出现逆温现象, 具有夜间稳定大气边界层温度分层效应。注意到相比14:00地面温度从288 K降到282.5 K, 降幅有5.5 K, 但是在高处温度仅由286 K下降到285 K, 仅下降了1 K。地表降温使得气流过山后减速时不再增温, 垂直方向上的逆温也抑制着垂直对流。高空处气温从14:00到24:00降幅很低是由于大气从白天的不稳定状态变为稳定状态, 而高空混合层不易受地面降温的影响, 在夜间转为了残留层, 致使其日变化值很小。
图9 24:00风速风向(a)、 气温(b)的侧边界条件与算例2(c)(d), 算例1(e)的模拟结果

Fig.9 The lateral boundary conditions of wind (a) and air temperature (b) and the simulation results of example 2 (c) (d), and 1 (e) at 24:00

图9(e)显示的是域内算例1模拟24:00的风速。和图9(c)相比, 空间上流场和风速的分布更为平直, 气流在过山后略有抬升, 和图9(c)相比并没有太多不同之处。
图10(a)、 (b)是算例2模拟14:00和24:00地面及XZ截面的湍动能分布。从图10(a)、 (b)中可以看出, 14:00湍动能在空间上的分布很不均匀, 最大值出现在山顶上空和山后山空, 24:00垂直方向上湍动能远小于14:00。白天垂直方向上分布的湍动能远大于夜间, 这说明大气湍流在白天不稳定状态下受浮力作用而增强, 夜间因为太阳辐射消失, 浮力作用减弱乃至消失使得湍流减弱, 夜间垂直方向上出现的逆温也是大气湍流减弱的原因之一。
图10 算例2、 算例1模拟的湍动能分布

Fig.10 The simulation result of example 2 and 1 for turbulent kinetic energy distribution

图10(c)、 (d)是算例1模拟14:00和24:00地面及XZ截面的湍动能分布。图10(c)、 (d)与算例2模拟的结果相比湍动能要小很多, 尤其是在14:00, 图10(c)没有模拟出白天不稳定状态下浮力效应对大气湍流的增强。算例1模拟的14:00和24:00的湍动能分布没有显著的区别。从两图中可看出, 在垂直方向上山地区域和山后附近区域上空的湍动能比周围区域要大, 仅能表明气流在过山时因受到山体结构的阻挡, 在受迫抬升后湍流强度有所增强。

4 结论与讨论

在Fluent模式的基础之上修改传统的RANS湍流参数化方案, 选用大涡模型（LES）, 并通过Fluent的UDF接口引入温度层结的浮力模型, 得到修改后的Fluent模式。将WRF与Fluent模式耦合, 针对吉山测风塔地区构建出两个算例, 其中算例1中Fluent仍使用传统RANS湍流方案（realizable k-ε 方案）, 不引入外部浮力模型, 算例2使用修改后的Fluent模式。使用算例1、 2和WRF一起模拟吉山站附近地形的流场和气温场。
通过使用观测资料对比2010年12月至2011年3月4个月吉山测风站风速风向算例1、 2的模拟结果, 可知: WRF在耦合Fluent之后, 对吉山单站点风速风向的模拟效果得到改善。相比算例1, 算例2可以更精确地模拟吉山站各高度处的风速风向。2010年12月至2011年3月4个月70 m高度平均风速WRF、 算例1和算例2的平均误差分别为3.963, 2.727和2.224 m·s-1, 6层高度(4, 10, 30, 50, 70和100 m)平均风速算例1的相对误差分别为30.9%, 24.7%, 14.7%, 16.5%, 16.3%和10.2%, 算例2的相对误差分别为26.9%, 22.5%, 12.8%, 12.3%, 13.3%和9.3%。算例2与观测的70 m高度风向玫瑰图匹配程度总体上比算例1更高。
通过对比分析2011年3月6日14:00和24:00两个时刻算例1、 2两者模拟的流场和气温场, 可知: 算例2可以模拟出白天不稳定状态下的热浮力效应（比如气流过山后减速增温形成的热对流）和空间分布不均匀且较大的湍动能, 能模拟出夜间大气边界层的逆温现象, 而算例1则不能。以上结论体现了算例2的优越性。
对于WRF耦合Fluent两个算例模拟结果的验证, 需要更多站点和更高采样频率的观测资料加以验证, 本文受条件所限, 只使用了鄱阳湖地区单个站点的平均风速资料。而Fluent在应用于对近地层大气流场的模拟时, 本身也存在一些没有解决的问题, 比如Fluent设置中下垫面参数过少以及水汽过程的耦合模拟问题等等, 这些都是将来进一步工作的方向。
综上所述, 本文修改了现有的Fluent模式, 改善了WRF耦合Fluent模式的模拟方法对江西吉山测风塔附近近地层风场的模拟效果, 为大气边界层底层大气流场预测和多尺度耦合模拟提供新的方法和依据, 也为未来发展和改进Fluent模式打下了基础。

感谢王咏薇老师和程雪玲老师对我的精心指导, 感谢鄱阳湖吉山测风塔观测人员的数据采集工作。

参考文献

原文顺序 | 文献年度倒序 | 文中引用次数倒序
Bechmann A, Berg J, Courtney M S, al et, 2009.The Bolund experiment: Overview and background.[R].Risø-R-1658(EN).

Beljaars A C M, Walmsley J L, Taylor P A, 1987.A mixed spectral finite difference model for neutrally stratified boundary-layer flow over roughness change and topography[J].Boundary-Layer Meteorology, 38(3): 273-303

Benchmann A, 2006.Large eddy simulation of atmospheric flow over complex terrain[D].Technical University of Denmark, 1-107.

Chow F, Street R, 2009.Evaluation of turbulence closure models for large-eddy simulation over complex terrain: Flow over Askervein Hill[J].Journal of Applied Meteorology & Climatology, 48 (5): 1050-1065.

Coppin P, Bradley E F, Finnigan J J, 1994.Measurements of flow over an elongated ridge and its thermal stability dependence: The mean field[J].Boundary-Layer Meteorology, 69(1/2): 173-199.

Hunt J, Tampier F, Weng W S, al et, 1991.Air flow and turbulence over complex terrain: A colloquium and a computational workshop[J].Journal of Fluid Mechanics, 227(6): 667-688.

Kim H G, Patel P C, Lee C M, 2000.Numerical simulation of wind flow over hilly terrain[J].Journal of Wind Engineering and Industrial Aerodynamics, 87: 45-60.

Koblitz T, Bechmann A, Sogachev A, al et, 2011.Modification of CFD code to model the non-neutral atmospheric boundary layer[R].13th International Conference on Wind Engineering,Amsterdam,Netherlands.

Kristof G, Racz N, Balogh M, 2007.Application of anasys-fluent on mesoscale atmospheric flow simulations[R].ANSYS Conference & 25th CADFEM User’s Meeting,November,Congress Center Dresden,Germany

Launder B E, Spalding D B, 1972.Lectures in mathematical models of turbulence[M].London: Academic Press, 1-7169.

Leblebici E, Ahmety G, Tuncer I H, 2011.Terrain modeling and atmospheric turbulent flow solutions based on meteorological weather forecast data[R].6th Ankara International Aerospace Conference,METU,Turkey Ankara

Li L, Hu F, Cheng X L, al et, 2006.Numerical simulation of the flow within and over an intersection model with Reynolds-averaged Navier-Stokes method[J].Chinese Physics B, 15(1): 149-155.

Meissner C, Schmitt C, Goretzki B, 2011.Non-neutral wind conditions in complex terrain.[J].Geophysical Research Abstracts, 13: EGU2011-12722.

Mortensen N G, Landberg L, Troen I, al et, 1993.Wind Altas Analysis and Application Program (WAsP): Vol.2: User's guide[M].Risø National Laboratory, 1-133.

Russell A, 2009.Computational fluid dynamics modeling of atmospheric flow applied to wind energy research[D].Boise State University, 1-85.

Shih T H, Liou W W, Shabbir A, al et, 1995.A new k-ε eddy viscosity model for high Reynolds number turbulent flows: Model development and validation[J].Computers & Fluids,24(3): 227 -238.

Shi R F, Cui G X, Wang Z S, al et， 2008.Large eddy simulation of wind field and plume dispersion in building array[J].Atmospheric Environment, 42(6): 1083-1097.

Stangroom P, 2004.CFD modeling of wind flow over terrain[D].The University of Nottingham, 1-299.

Sumner J, Watters C S, Masson C, 2010.CFD in wind energy: The virtual, multiscale wind tunnel[J].Energies, 3(5): 989-1013.

Taylor P A, Teunissen H W, 1987.The Askervein Hill Project: Overview and background data[J].Boundary Layer Meteorology, 39(1/2): 15-39.

Thunis P, Bornstein R, 1996.Hierarchy of mesoscale flow assumptions and equations[J].Journal of the Atmospheric Sciences, 53(3): 380-397.

Tseng Y H, Meneveau C, Parlange M B, 2006.Modeling flow around bluff bodies and predicting urban dispersion using large eddy simulation[J].Environmental Science & Technology, 40(8): 2653-2662.

U.S, National Renewable Energy Laboratory, 2006.Technical Report-China wind energy resource assessment[M/OL].http: //swera.unep.net.

Wang W G, Shaw W J, Seiple T E, al et, 2008.An Evaluation of a Diagnostic Wind Model (CALMET)[J].Journal of Applied Meteorology and Climatology, 47(6): 1739-1756.

Weng W S, 1997.Stably stratified boundary-layer flow over low hills: A comparison of model results and field data[J].Boundary-Layer Meteorology, 85 (2): 223-241.

Wilcox D C, 1998.Turbulence modeling for CFD[M].DCW Industries, 1-543.

Wood N, 2000.Wind flow over complex terrain: A historical perspective and the prospect for large-eddy modelling[J].Boundary-Layer Meteorology, 96 (1/2): 11-32.

Yakhot V, Orszag S A, 1986.Renormalization group analysis of turbulence[J].Journal of Scientific Computing, 1(1): 3-51.

Yu W, Benoit R, Girard C, al et, 2006.Wind Energy Simulation Toolkit (WEST): A wind mapping system for use by the wind energy industry[J].Wind Engineering, 4034(1): 15-33.

曹杨, 陈洪滨, 王普才, 2017.声雷达资料可靠性及近地面边界层风场特征分析[J].高原气象, 36(5): 1315-1324.DOI: 10. 7522/j.issn.1000-0534.2016.00100.

程雪玲, 胡非, 2005.大气边界层内羽流扩散研究[J].力学学报, 37(2): 148-156.

程雪玲, 胡非, 2006.复杂地形网格生成研究[J].计算力学学报, 23(3): 313- 316.

方艳莹, 徐海明, 朱蓉, 等, 2012.基于WRF和CFD 软件结合的风能资源数值模拟试验研究[J].气象, 38(11): 1378-1389.

何晓凤, 周荣卫, 朱蓉, 2010.MM5与CFD软件相结合对复杂地形风资源模拟初探[J].资源科学, 32(4): 650-655.

姜平, 刘晓冉, 朱浩楠, 等, 2019.复杂地形下局地山谷风环流的理想数值模拟[J].高原气象, 38(6): 1272-1282.DOI: 10.7522/j.issn.1000-0534.2019.00019.

李磊, 张立杰, 张宁, 等, 2010.FLUENT在复杂地形风场精细模拟中的应用研究[J].高原气象, 29(3): 621-628.

李晓霞, 黄涛, 王兴, 等, 2017.兰州新区近地层风场时空特征分析[J].高原气象, 36(4): 1001-1009.DOI: 10.7522/j.issn.1000-0534.2016.00092

刘丽珺, 梁友嘉, 2018.集成CFD与Kalman滤波的微尺度风电场风功率预报方法[J].高原气象, 37(4): 1061-1073.DOI: 10. 7522/j.issn.1000-0534.2017.00098.

刘郁珏, 苗世光, 胡非, 等, 2018.冬奥会小海坨山赛区边界层风场[J].高原气象, 37(5): 1061-1073.DOI: 10.7522/j.issn.1000-0534. 2018.00034.

肖仪清, 李朝, 欧进萍, 等, 2009.复杂地形风能评估的CFD方法[J].华南理工大学学报（自然科学版）, 37(9): 30-35.

周荣卫, 何晓凤, 2018.新疆哈密复杂地形风场的数值模拟及特征分析[J].高原气象, 37(5): 1413-1427.DOI: 10.7522/j.issn. 1000-0534.2018.00021.

Options
文章导航

/

〈 〉

本系统由北京玛格泰克科技发展有限公司设计开发 技术支持：support@magtech.com.cn