为了有效应对气候、环境变化和能源短缺，最近20年来世界上风能资源的开发利用得到了迅猛发展。我国幅员辽阔，风能资源极为丰富。随着风电产业的发展，对风能资源评估的精度要求越来越高，一般都采用数值方法分析风资源分布。国内外对复杂地形风场结构的数值模拟有大量的研究。归纳起来可分为诊断模式(diagnostic model)和预报模式(prognostic model)两大类。诊断模式是从实验观测数据出发构建稳定状态的流场，方程中不含时间项；预报模式是从初始流场开始，同时考虑边界条件，得到气象场的时间演化，方程中含有时间项。值得注意的是，如果预报模式用来计算定常过程，可以当作诊断模式。

诊断模式又可分为两种，一种是基于线性化方法简化Navier-Stokes方程的定常解，称为线性模 式(linearized model)；另一种是根据质量守恒定律分析气象场数据，称为质量守恒风场调整模式(mass consistent model)。线性模式最早是Jackson and Hunt(1975)提出来的，他们采用线性化和混合长理论分析了与二维小山的大小和形状有关的越山气流，将边界层分为4个子层，给出了相应的速度分布公式，但该结论有效的前提是需要山体坡度很小。随后，Hunt et al.(1988)进一步分析了剪切流越过二维和三维小山的情况，并放宽了来流风速剖面和粗糙度的变化要求，但研究对象仍然仅限于缓坡度的小山，即要避免分离的产生。20世纪80～90年代，丹麦Riso国家实验室在其理论基础上，发展了一个用于风电场微观选址的资源分析工具软件——WasP(Wind Atlas Analysis and Application Program)，该软件适合较为平坦的地形、小范围的风资源调查(Mortensen and Landberg，1993)。Kaimal and Finnigan(1994)研究认为此线性理论在计算未发生流动分离时的最大加速效应误差可以小于15%～18%。Beljaars et al.(1987)也在此线性理论的基础上开发了边界层数值模拟程序MSFD模式(Mixed Spectral Finite Difference model)，Taylor(1998)、Weng et al.(2000)先后运用其研究了山体几何尺寸和地面粗糙度对加速比的影响，提出了可以考虑地面粗糙度影响的山顶加速比算法，称之为“新算法”(new guidelines)。质量守恒风场调整模式通过近地面和探空观测资料构建三维风场。首先将观测资料插值到计算格点上，然后调整插值结果以满足质量守恒。其中，插值方法非常关键，相同的观测资料由于插值方法的不同，最后得到的风场会相差很大。而对插值结果的调整则要根据研究问题的不同而进行，比如对局地风场，可根据不可压连续方程调整；对污染物长程输运问题，则采用可压连续方程进行调整(Ishikawa，1994)。美国环保局推荐的空气质量扩散模式CALPUFF(California Puff Model)中的气象模块CALMET(California Meteorological Model)就是利用质量守恒原理对风场进行诊断的模式，包含了客观化的参数分析、陡坡地形的参数化处理、地形影响下的动力学流体效应、特殊地形对大气流体的阻滞效应、辐散散度最小化处理，以及专门为处理海陆边界层和大面积水体区域上空的气体扩散的微气象学处理算法(Wang et al.,2008)。

预报模式则是通过质量、动量、热量、水汽以及其他如气溶胶等守恒方程的时间、空间积分预报大气系统的演化过程。各种预报模式之间最主要的区别在于守恒方程采用的近似方案，比如不可压、Boussinesq假设、静力平衡等，以及对物理过程的参数化方案，比如云过程、沉降过程、湍流、通量等。各近似方案和参数化方案都有自己的局限性(Thunis and Bornstein，1996)。此外，不同模式数值求解非线性偏微分方程的方案也多种多样，包括网格划分、数值方法、初边界条件设置、坐标系的选择等均不相同。目前被广泛使用的预报模式有WRF、MM5、RAMS、ARPS等，其近地层、边界层、次网格参数化方案均对近地层风速模拟结果有较大影响，对于复杂地形则更多地采用的是地形追随、sigma和z*坐标系，而非笛卡尔坐标系，通过方程组的坐标变换来描述复杂地形，在数值计算方法上以差分格式为主，且在模拟的前处理阶段需要对地形进行不同程度的平滑，以获得计算上的稳定性，在近地面获得更高的精度，当遇到极为陡峭的地形时，有可能出现积分溢出的情况(Finardi et al.,1997)。

为了满足复杂地形风场的精细化数值模拟，单纯用预报模式很难达到需求，即使能够实现精细化模拟的预报模式，对于区域性模拟需耗费大量机时，因此，需要通过动力降尺度的方法来实现精细化数值模拟，出现了预报模式和诊断模式相结合的办法计算三维流场。美国True Wind Solutions公司将中尺度数值模式MASS(Mesoscale Atmospheric Simulation System)与小尺度质量守恒的风场线性模式WindMap结合，发展成MesoMap应用于风能资源评估(Brower et al.,2001)，分辨率达到100～1000 m。加拿大气象局将中尺度模式MC2与小尺度模式Ms-micro相结合建立了精细化风能资源数值模拟评估模型WEST(Wind Energy Simulation Toolkit)，并绘制了5 km×5 km分辨率的加拿大全国风能资源图谱(Yu et al.,2006)，张德等(2008)用该模式对我国风能资源进行了模拟，其模拟结果与气象台站观测资料对比分析的结果表明，WEST在中国中东部海拔较低地区的模拟能力较好，但对于中国西部高原地区的模拟结果不理想。周荣卫等(2010)利用中尺度气象模式MM5和微尺度模块Calmet组合建立了一个适用于区域性风能资源评估的数值模式系统，用该模式系统分别对代表我国沿海风能丰富区的江苏省沿海地区和复杂地形条件下的北部风能丰富区的甘肃省酒泉地区的风能资源进行1 km×1 km 高分辨率的逐时模拟，并通过与模拟范围内测风塔观测资料的对比检验该模式系统对我国不同地形条件下的风能资源评估能力，结果表明该模式系统对甘肃省酒泉地区的模拟值与实测值的对比结果要优于江苏省沿海地区的对比结果。

近年来，随着计算机能力的增强，以往用于空气动力学精细流场计算的CFD模式越来越多地在气象领域得到应用，比如用于计算城市微尺度风场及污染扩散(Cheng and Hu，2005；Neophytou and Britte，2005)。由于具有成熟的网格生成模块，能够自适应地生成各种复杂地形上的贴体网格，以及强大的后处理模块显示局部的复杂流动，CFD模式更适应于复杂地形条件下的边界层精细流场的模拟。程雪玲和胡非(2005，2006)利用地形高程的格点数据建立地形的计算网格模型，输入CFD模式模拟大气边界层内复杂地形的羽流扩散。Stangroom(2004)利用CFX-5对Askervein山的风场进行数值模拟，结果优于WAsP，与实验值吻合很好。Bechmann(2006)基于EllipSys3D发展了大涡模拟(Large Eddy Simulation，LES)模式，同样模拟Askervein山的绕流问题，得到的平均风速和湍流脉动与实验均吻合。李磊等(2010)用商业CFD软件Fluent 模拟了复杂地形的精细风场。在此基础上，人们开始研究用中尺度预报模式和CFD模式结合进行复杂地形风场的数值模拟(何晓凤等，2010；方艳莹等，2012；Li and Chan，2012；Li et al.,2013)。

中尺度大气模式通常采用地形追随坐标，通过方程组的坐标变换来描述复杂地形，在数值计算方法上以差分格式为主，且在模拟的前处理阶段需要对地形进行不同程度的平滑，以获得计算上的稳定性(Dudhia et al.,2005；Walko and Tremback，2006)，当遇到极为陡峭的地形时，有可能出现积分溢出的情况，无法解决复杂地形的风能资源评估。

基于CFD方法的数值模式一直是用来进行空气动力学研究的，更适应于复杂地形条件下的边界层精细流场的模拟，2011年Riso实验室组织了一次实验，在长130 m、宽75 m的Bolund岛上，安装了11台测风仪，实验进行了30 d。这次实验对不同机构57个CFD模式，包括线性类CFD，如WAsP、WindFarmer、雷诺平均Navier-Stokes方程(RANS)类CFD、LES类CFD等进行了盲比，对平均风速的计算结果比较后，发现线性CFD模式差，RANS类CFD模式最好(Bechmann et al.,2011；Berg et al.,2011)。Fluent是国际知名的CFD模式，集成了非常丰富的参数化方案和物理模型，以及强大的建模和网格划分功能，适用于模拟复杂地形的风场。Cabezón et al.(2007)用Fluent和WAsP计算了高度复杂地形的风场，前者的计算精度更高。Sadek et al.(2011)用Fluent和线性模式Flowstar计算了复杂地形的风场和浓度扩散，Fluent的计算结果更准确，特别是对于起伏的山地和陡坡形成的回流区。Mandas et al.(2004)用Fluent以Askervein山为例对其计算精度进行检验，其结果要优于其他模式。

利用Fluent模拟复杂地形的风场，在边界条件的处理上缺乏灵活性，通常只能假设边界上入流条件保持幂指数风廓线状态，无法保证数值模拟中 能正确地包含真实的中尺度环流信息。Li et al.(2004)用RAMS+Fluent构造耦合模式系统成功地用于城市冠层气象场方面精细模拟的模拟研究。本文的耦合模式系统采用WRF模式，通过嵌套网格到某一最小内层尺度(一般是几公里)，然后通过耦合CFD模式Fluent获得高水平分辨率(30～100 m)和高垂直分辨率(150 m以下每10 m一层)的风速分布资料。得到的资料还将与观测站点(气象塔)的实测资料进行对比，检验该系统的误差。

在计算风电场尺度的风场时，Fluent假设空气运动是三维不可压缩流，忽略地转科里奥力的影响，且不考虑热对流的影响，所有湍流都是因为机械作用产生的(Versteeg and Malalasekera，1995)。

空气运动可以分解为均值u_i和脉动u'_i：

同样气压也这样分解：

Navier-Stokes方程可以重新写为

方程(3)和(4)就是RANS方程，其中，t为时间，u_i、u_j分别为i、j方向速度，d_ij为Kronecker delta函数，ρ为空气密度，u为动力粘性系数。方程(4)中多了一个未知项——雷诺应力张量，代表速度脉动量的相关性，需要被参数化。

有几种方法来参数化雷诺应力，如k−w、k−e、realizable k−e方案、RNG k−e方案、雷诺应力输送模型以及LES模型。LES模型能够给出最好的参数化结果，但太耗计算资源。Fluent 6.3集成了3种基于Reynolds平均思想的湍流雷诺应力方案：标准的k−e方案(Launder and Spalding，1972)、RNG k−e方案(Yakhot and Orszag，1986)和realizable k−e方案(Shih et al.,1995)，其中realizable k−e湍流闭合方案被广泛地应用在风工程、污染物扩散和城市微尺度环境的研究，此方案适当地调和计算精度和计算资源，确保可靠的结果和稳定。本文采用中尺度气象模式耦合CFD软件Fluent的方法研究复杂地形的风场。

当前，国际上关于复杂地形风场的模拟(或预报)大都采用的是中尺度气象模式(例如MM5、WRF等)，通过嵌套网格到某一最小内层尺度(一般是几公里)，然后通过动力降尺度(例如使用CALMET模式)或统计降尺度的方法，以及大气边界层风廓线分布理论或半经验公式，获得高水平分辨率(几百米到1 km)和高垂直分辨率(一般要求150 m以下每10 m一层)的风速分布资料。这样得到的资料要与观测站点(气象塔)的实测资料进行对比分析，从而改进模式和模拟方法，最终得到满足一定精度要求的数值模拟结果和风能资源分布图。本文则是把上述国际上流行方案中的动力降尺度部分，改为Flunet精细模拟。

鄱阳湖平原位于江西省北部、平原外侧，低丘岗地广布，地面波状起伏，海拔在50～100 m 之间，内侧湖滨圩区，海拔多在 20 m 以下，地势低平，港汊纵横草洲滩地连片，鄱阳湖位居中央(见图 1)。从湖口到吉山形成90 km长的狭管湖道，南部大湖体水面广阔，湖区周围地形复杂，北侧紧邻庐山，东侧为江西最高峰黄岗山海拔2157 m，其东北、西北和西侧均为山地。由地形、山谷风、湖陆风的相互作用，形成复杂的大气边界层结构、局地气候和风能资源的分布。在其湖道地区，由于风向稳定、风力资源优质，建有多个风电场和测风塔。

图 1

Fig.1

图 1 鄱阳湖地区地形 Fig.1 Terrain of the Poyang Lake area

首先使用WRF模式进行6 km分辨率的模 拟；然后将WRF模拟的结果作为Fluent的边界条件，进行CFD计算得到测风塔周围6 km×6 km范围内水平50 m分辨率、150 m高度以下垂直10 m分辨率的风能资源分布结果。WRF采用3层网格嵌套，网格格距分别为54 km、18 km、6 km(如图 2所示)。

图 2

Fig.2

图 2 WRF网格划分 Fig. 2 Grids of WRF model

鄱阳湖地区天气系统复杂，对风场产生影响的环流系统较多：冬季常受西伯利亚(或蒙古)冷高压影响，盛行偏北风，春季受季风影响，夏季多为副热带高压的控制，盛行偏南风(见图 3)，有时受台风影响，且局地强对流天气时常发生。图 4给出了沙岭观测塔70 m高度上从冬季西北大风天气到夏季偏南风时模拟风速和观测风速的对比，可见两者变化是一致的。图 5～8给出了其他4座测风塔各层高度上从冬季到夏季模拟与观测风速月平均误差，可以看到冬季由于大风天气较多，经常会出现寒潮、冷峰等转折型天气，模拟的误差较大，而夏季风速平稳，模拟的结果要相对较好，这也能够满足风电场风功率预报的要求，大风天气由于超过风电曲线的截止风速，准确预报大风并不是必须的，而只有准确预报截止风速以下的风速，才能保证风电预报的准确性。

图 3

Fig.3

图 3 微尺度风场的水平和垂直网格 Fig. 3 Horizontal and vertical grids of the flow field

图 4

Fig.4

图 4 沙岭测风塔70 m高度上模拟与实测风速比较Fig. 4 Comparison of simulated and observed velocities at 70-m height at Shaling tower

图 5

Fig.5

图 5 矶山月平均风速模拟与实测误差Fig. 5 Errors of monthly average velocity between simulation and observation at Jishan tower

图 6

Fig.6

图 6 同图 5，但为吉山Fig. 6 Same as Fig. 5，but for Jishan tower

图 7

Fig.7

图 7 同图 5，但为皂湖Fig. 7 Same as Fig. 5，but for Zaohu tower

图 8

Fig.8

图 8 同图 5，但为狮子山Fig. 8 Same as Fig. 5，but for Shizishan tower

鄱阳湖地区地形复杂，地面条件如地形起 伏、水陆交界、地表粗糙度或者降水和云量的不均匀分布，都可能造成中、小尺度大气运动状况的改变，这种情况下的风场模拟，地形的影响比较关键。WRF模式模拟所采用的地形地表资料是USGS24 类 30 s分辨率资料，格点间隔为 6 km，Fluent所采用的地形资料是GDEM 水平分辨率30 m资料，格点间隔为50 m。WRF模拟的地表风速往往低估山峰处的风速，而高估平原处的风速。图 9和图 10是耦合模式得到的离地面70 m高度上的风速的矢量图和等值线图，从图上可以看到，山上的风速要远高于山下，且从矢量图上可以很明显地看到风翻越山顶时的风向变化，得到了较为合理的当地风场的分布，对风能资源普查有参考价值。

图 9

Fig.9

图 9 70 m高度上流场风速矢量图Fig. 9 Wind speed vectors at 70-m height

图 10

Fig.10

图 10 同图 9，但为风速等值面Fig. 10 Same as Fig. 9，but for wind speed contours

云南省陆良县板桥镇杨梅山风电场位于(25°6.95′N，103°52.171′E)，平均海拔2200 m。图 11a中为风电场中66 台风机的位置，图 11b是根据当地地形建立的精细网格，其中的绿线为风机所在的位置，红线是测风塔所在位置。

图 11

Fig.11

图 11 杨梅山风电场(a)实际地形和(b)精细网格Fig. 11(a)Complex terrain and (b)fine grid representation at Yangmeishan wind farm

图 12是66台风机轮毂高度的风速模拟结 果，图 13是模拟测风塔上的风速与观测结果的对比，并与相似Kalman滤波(MOS)订正的结果 进行了对比，相似Kalman滤波适用于天气系统 剧烈转变的情形，可以看到数值模拟的结果、订正结果与观测的变化趋势是一致的，但在风速发生突变的时候两者均有误差。因此，可考虑将该耦合 系统与统计预报模式组合，在进行风电预报时针对两种预报方法在不同的天气过程和不同风电场区域的预报表现，进行细致分类，建立集成预报系 统。实际预报时根据分类选择统计预报还是动力预报。

图 12

Fig.12

图 12 66座风机轮毂高度的风速Fig. 12 Wind speeds of the 66 turbines at the hub height

图 13

Fig.13

图 13 风速预报结果比较Fig. 13 Comparison of wind speed between different forecasting results

本文通过中尺度气象模式和计算流体力学模式耦合进行复杂地形风场精细模拟，得到水平分辨率50 m、150 m以下垂直分辨率10 m的精细风场，风机轮毂高度上的风速变化，同时比较了测风塔上模拟与观测的结果，进行了月平均的误差分析，效果令人满意。这两个模式是完全不同两个独立的模式，各自有各自的优点。风场范围内的大气运动不仅受区域中尺度环流的影响，还与局地小尺度运动有关。耦合这两个模式真正考虑到了大尺度因素(全球模式提供边界条件给中尺度模式)、中尺度因素和小尺度因素对某个具体局地环境的影响，有巨大的应用前景，例如风能资源评估、城市污染预报、核扩散评估和大型建筑群小环境设计等。