基于自主软件LSWS 的典型地形大气CFD 模式研究与应用

李树民，孙 壮，李泽祥，李鑫洋，罗朝匀，李 辉，孙鑫宇，刘伟毅

(成都流体动力创新中心，成都 610072)

0 引言

风能资源作为最具开发潜力的可再生能源之一，其精细化评估直接关系到风能利用效率和极端风速下的风电机组安全性。近年来，风电机组不断向大型化方向发展，叶尖高度超过200 m、风轮直径超过150 m 的大型风电机组屡见不鲜，这种发展趋势给风能资源评估带来了新的挑战。现有的风资源评估技术建立在经典的近地层相似理论基础上，仅适用于高度100 m 左右的近地层[1]。在100～300 m 高度的大气边界层（atmospheric boundary layer,ABL）下层不再是地表摩擦力和气压梯度力二力平衡的大气运动，而是加入了地转偏向力的三力平衡情况。因此重新认识风特性，拓展或重构风资源理论和方法，揭示典型地形影响下300 m 高度内的风和湍流特性及其形成机理，已成为国内外风电领域的新的研究重点。计算流体力学（CFD）方法可以获得精细化的流场细微结构，已成功应用于很多工业领域，在精细化风资源评估中具有一定潜力。由于目前传统工业用CFD 软件尚不具备大气效应的求解能力，而主流的风资源评估软件均被国外垄断，使得我国风资源精细化评估工作受到严重制约，因此建立一套能够适用于我国典型地形风场特征的大气模式的仿真模拟方法，使其具备中性、非中性（稳定、不稳定等）、地转效应等大气风场仿真能力，将中尺度大气模式模拟方法移植到传统CFD 中（本文中称之为大气CFD 模式），实现天气背景风场的大气模式与传统CFD 模拟的有效衔接，将对我国风资源精细化评估、微观选址、风电预测起到重要支撑作用。

风资源评估经历了从无热力变化的中性大气到考虑热力变化的非中性大气的研究历程。在中性大气方面，热力作用不明显，风速较为强劲，研究者多集中在计算风工程领域。Richards等[2]最早提出了中性大气边界层流场模拟方法及边界条件，并在计算风工程领域得到广泛应用。Hargreaves等[3]通过改变壁面边界条件，并在计算域上部考虑大气的剪切作用，实现了中性风速廓线在较长计算域中的保持性研究。杨易等[4]从边界条件应满足湍流模型的角度出发，引入了新的入流边界条件，在原有速度入口廓线的基础上，给出了入口湍动能边界条件随高度变化关系式。张晓东[5]针对壁面函数模型与入流剖面不一致会导致流向上出现不合理的梯度变化情况，提出了一种壁面函数模型，缓解了入口速度廓线在计算域中的梯度变化，并应用于Bolund Hill 和Askervein Hill两个典型山地风场的仿真模拟与验证。遆子龙等[6]基于管道流动的粗糙度理论研究了地表粗糙度对山区峡谷桥梁区域风场的影响。

在非中性大气风场模拟方面，研究主要集中在风能领域，以丹麦科技大学（DTU）的研究为代表[7]，主要根据莫宁-奥布霍夫相似理论，将莫宁-奥布霍夫长度、热浮力、科氏力等大气特征参数引入湍流模型，实现了非中性大气的仿真模拟，在这一思想指导下开发出了EllipSys3D软件，并衍生出了风资源评估软件WAsP[8]。法国美迪公司的Meteodyn WT 风资源评估软件采用了同样的思想，基于混合长度理论将莫宁-奥布霍夫长度引入一方程湍流模型中，并通过入口速度和温度廓线的分层理论，实现了非中性大气的仿真[9]。在复杂地形非中性大气仿真方面，比较有代表性的是采用域前模拟方法来开展[7]：第一步，在地形计算域前建立一个空的计算域，利用循环计算的方法，配合非中性的湍流模型和边界条件，发展出平坦地形条件下的非中性风场状态；第二步，以第一步发展出的风速为边界条件，施加进带有复杂地形的计算域中，实现非中性风场的仿真模拟。虽然经过了多年的发展，但采用CFD 方法模拟大气风场依然是个非常有挑战性的难题[10-11]。

由于国外商业软件的通用性、便捷性和强大的市场运作能力，严重压缩了我国自主软件的生存空间。目前在风电领域，国内已有的风电场数值模拟软件主要依靠国外。国产自主软件大多是为航空、航天、航海等问题计算而开发，对10 km 级及以上的大范围风场仿真的工程应用还比较少见。我国依托数值计算开展的风力发电技术创新研究也面临着CFD 软件使用受限的风险。鉴于此，本文避开通用软件市场，在专业软件方面精准发力，开发了能够模拟仿真我国不同地形、不同大气稳定度特征的自主软件LSWS，并开展了典型地形的仿真验证工作。本文所述典型地形是指我国具有较强风资源开发潜力区域所包含的主要地形，包括：内蒙古草原、丘陵、长江中下游平原、山地、秦岭山脉、高原山地等。在研发大气CFD模式的过程中，选取了典型地形中的两种地形开展研究，包括内蒙草原锡林浩特平坦草原地形和山西丘陵地形。在这两种地形下，我国建有大量的风电场，对于研究风特性、风能特性有重要的参考意义。文章第一节为大气CFD 模式及边界条件，详细介绍了大气CFD 模式理论，基于锡林浩特风场实测数据提出的修正Gryning 风廓线模型和拟合廓线模型，以及所用计算方法和自主软件架构；第二节为我国典型锡林浩特平坦草原地形、山西丘陵地形风场计算结果，并与商业软件Fluent 的计算结果进行了对比；第三节给出了结论及建议。

1 理论及方法

1.1 大气CFD 模式及边界条件

1.1.1 大气CFD 模式

为了实现对大气边界层风场的模拟，考虑不同大气稳定度的影响，在质量守恒方程和动量守恒方程的基础上引入了位温方程。控制方程组可以写为：

其中：ui为i方向的平均速度分量；µ为空气动力黏度；µt为湍流涡黏；P为压力；ρ 为密度；σθ为位势温度θ的施密特数，θ的表达式为 θ=T(P0/P)R/cp；T为大气温度；P0≈1×105Pa 为标准大气压；R=8.313 J/(mol·K)为普适气体常数；cp≈1 J/(g·K)为干空气在恒定压力下的比热容。

大气边界层内部的大尺度运动受到热分层和科里奥利效应的影响，分别由地球自转和惯性质量引起。当对整个ABL 进行建模时，这些效应可以通过附加的源项（Sv）引入RANS 方程组：

其中：gi为重力加速度，gi=(0,0,-g)T；ρ0为参考密度；εi=(-1,1,0)T；fc=2ΩEsinλ为科里奥力参数，ΩE为地球自转速度，λ为纬度。

1.1.1.1 中性大气湍流闭合方程

使用k-ε湍流模型封闭上述控制方程组。适用于中性ABL 的k-ε模型表达式为：

其中：σk和 σε分别为k和 ε 的施密特数；Cε1和Cε2为模型系数；Pk为湍动能的生成项。通过求解湍动能k和耗散率 ε的两个偏微分输运方程，得到湍流涡黏。由此得到的混合长度lt和涡流黏度 µt用k和 ε可表示为：

其中Cµ为常数。

1.1.1.2 非中性大气湍流闭合方程

在非中性条件下，为了考虑热稳定性的影响，本文参考Koblitz[7]论文中的理论公式，对k-ε模型进行修正，最大的特色是通过经验的湍动能理论估计对全局最大混合长度进行显式积分计算，避免了原文中繁琐且没有理论依据的预前模拟才能求解的限制。

湍动能方程中增加的浮力源项B与系数Cε3一起出现在耗散方程。浮力源项B的表达式为：

在热不稳定条件下，B为正，支持湍动能的产生；而在稳定条件下，B为负，抑制湍动能的产生。

附加扩散项D的表达式为：

其中，系数Cµ为常数。

其中，le为全局最大混合长度。在中性条件下的平坦草原地形上，le=l0=0.000 27G/fc，G为地转风速。对于非中性条件和复杂地形，有：

通常取系数α=0.075，用于保证在中性条件下l0和lMY两/个长度尺度完全相同。对湍动能k取理论估计并忽略大气边界层高度以外的湍动能，这里u*为大气边界层的摩擦速度。

式（15）可近似表达为lMY=αhg/2 。hg=u*/(6fc)为大气边界层高度的估算。αB的表达式为：

其中，理查德森数Rig=B/Pk，B依赖于涡黏 µt、重力加速度gi以及密度 ρ。

在数值模拟过程中对多种湍动能入口廓线进行了对比。对湍动能廓线的一种拟合形式为：

其中，C1和C2为拟合系数。针对这一湍动能廓线形式，对lMY的积分表达式进行积分可得：

数值模拟系数取值见表1。

表1 数值模拟系数取值Table 1 Values of the numerical simulation coefficients

1.1.2 壁面函数

在地形表面附近，通过壁面函数来计算地面上方第一层网格中心高度的函数值。在本文计算中，热稳定度的影响通过入口廓线来体现，壁面函数的形式基于平衡假设来推导，唯一的湍流能量产生机制是剪切作用，不考虑浮力等因素的影响。对数速度剖面的表达形式为：

其中，E为对数分布律中的常数，通常取8.143。由式（19）可推导出相关物理量的壁面函数。速度的壁面函数通过壁面切应力 τ0来实现，

其中：kp为网格中心处的湍动能；up为网格中心处的速度；为地面到网格中心zp的无量纲距离，

第一层网格中心的温度通过壁面热流密度q0来计算，

其中：Tp为第一层网格中心处的温度；T0为壁面处的温度；λt为第一层网格中心与壁面之间区域的当量导热系数，且有关系式：

其中：σl为分子Prandtl 数（即Pr）；σT为湍流Prandtl 数。

1.2 锡林浩特风速测量实验

入口条件是影响大气CFD 模式准确性的重要因素。以内蒙古自治区锡林浩特大气边界层湍流观测实验场的实测风数据对大气CFD 模式的入口廓线进行标定。该实验场位于锡林浩特东北方向，距离锡林浩特约30 km，为高原草场，地形平坦，零星分布一些低山丘陵，可以代表内陆平坦地形。前期朱蓉等[12]已经对该地区的风能资源和气候特征做了详细研究。该地区设立了100 m 观测塔和70 m 观测塔各1 个，本文所用数据的观测时段为2009 年6 月至2010年7 月。测风塔外观如图1 所示。对100 m 高度以下的风速，采用测风塔实测风数据。对于100 m 以上高度的风速，测风塔难以测量，采用激光雷达测风仪进行测量。

图1 锡林浩特风速测量实验Fig.1 Xilinhot wind measurement experiment

1.3 入口边界条件

本节给出了中性、非中性风廓线的经典模型及修正方法，对比了两种经典模型描述我国典型地形风特性的优缺点，基于锡林浩特实测风数据对经典模型进行了修正，提出了一种修正模型，同时基于分段拟合方法提出了一种拟合廓线模型，这两种模型均可作为大气CFD 模式入口边界条件。

1.3.1 廓线修正方法

大气边界层可以垂直分为层流底层、地表层、埃克曼层3 个部分：最低部分称为层流底层，其高度等于气动粗糙度高度z0；z0以上是湍流充分发展的地表层，其垂直范围可能从20 m 到100 m 不等[5]，具体取决于空气的热分层，通常为整个大气边界的10%；地表层之上是埃克曼层，其高度可能超过1 000 m，具体取决于空气的稳定性、科里奥利参数和地表的粗糙度高度。

风速随着高度的垂直升高而升高，是最重要的大气特征。最经典的方法是使用对数函数描述。在地表层：

式（25）只适用于中性条件。对于非中性条件，可采用莫宁-奥布霍夫相似理论对式（25）进行修正，

对于不稳定条件，修正方程[13-14]为：

对于稳定条件，对数风廓线修正方程为[15-17]：

其中：a=5；A=1；B=2/3；C=5；D=0.35。

当从地表层到达艾克曼层时，地球自转产生的科氏力不可忽略，上述方法不能准确描述速度随高度的变化，需要对整个大气边界层的风廓线进行统一描述。目前常用的描述方法可分为两种。

1.3.2 EKM 模型

第一个方法是采用分段拟合，并确保在两个区域内能够光滑过渡[18]。该模型在下文中称为EKM 模型：

其中：ug为地转风速；zp为地表层厚度；α0为表面层风与地转风的夹角。对中性条件，

对非中性条件，

ψ(zp/L)在稳定和不稳定条件下不同。对于稳定条件，

式中，a的取值在4.7～5 范围内。对于不稳定条件，

1.3.3 修正Gryning 模型

第二个方法是修正混合长度对高度的依赖性，其由Gryning等[19]提出。他们用公式表达了与高度相关的混合长度，以限定其随着高度的生长，从而把对数法则式的有效性延伸到了表面层以上。对于稳定条件：

对中性和不稳定条件：

式中，zi为大气边界层高度。

其他几项可由以下公式获得：

式（36）、式（37）在下文中称为Gryning 模型。不同稳定度下，风速随高度变化的规律如图2 所示。图中方块为基于测风塔和激光雷达的实测数据，实线为EKM 模型和Gryning 模型的拟合曲线。可以看出，在100 m 高度以下时，两种模型的拟合曲线与实测值均能良好吻合，其中Gryning 模型稍优于EKM 模型。随着高度增加，尤其是在稳定条件下，拟合曲线与实测值之间的偏差增大。这是由于理论模型高估了速度随高度增长的趋势，导致其在更高的高度下的速度估值偏高。

图2 风速随高度变化规律Fig.2 Variation of wind speed with height

相比于EKM 模型，Gryning 模型结构更加简单。因此，基于实测数据，对方程（36）、方程（37）右侧的第二项和最后一项进行重新拟合修正，获得Gryning修正模型。对于稳定条件：

对中性和不稳定条件：

Gryning 修正模型与实测数据的对比如图3 所示，修正模型与实测数据吻合良好。修正模型较好地修正了原模型与实测值随高度增加而存在的偏差。因此，在保证Gryning 模型整体结构不变的情况下，仅需要对其相关常数项进行简单的修正，即可与实测值吻合。

图3 Gryning 修正模型与实测值的对比Fig.3 Comparison between the modified Gryning model and the measured value

1.3.4 拟合廓线模型

除使用理论模型外，基于测风塔和激光雷达的实测数据，结合分段函数拟合方法和理论公式，本文构建了不同热稳定度下的风速u随高度z变化的入口廓线模型，描述如下。

当热稳定度分级小于或等于3 时：

当热稳定度分级为4 或5 时：

其中：u*为通过实测获得的地表摩擦速度；u4为通过实验测量获得的4 m 处风速；z0为地表摩擦度；L为莫宁-奥布霍夫长度；κ=0.41为冯·卡门常数；ug为地转风速，

A、B与稳定度有关（Garratt，1994）：

c和c′为与热稳定度有关的常数，由表2 可获得。

表2 不同热稳定度下拟合参数Table 2 Fitting parameters under different thermal stability degrees

拟合曲线与实测数据对比如图4 所示，模型拟合曲线在不同稳定度下均能与实测数据良好吻合，该模型能够为典型地形非中性大气CFD 模式计算提供准确的速度入口条件。

图4 拟合模型曲线与实测数据对比Fig.4 Comparison between the fitting model curve and the measured data

1.3.5 入口湍动能

在大气边界层流动中，被广泛采用的湍流动能k和动能耗散率ε的廓线可由以下两个方程给出：

其中，Cµ一般取0.03。本文基于Zhang[5]的工作，引入一随高度变化的修正项，修正上述湍动能方程：

图5（a）给出了实测数据与方程（53）曲线的对比。由图可见，仅在稳定条件下，方程（53）的预测值与实测值吻合良好；在弱稳定条件和不稳定条件下，方程（53）低估了湍动能；在中性和弱不稳定条件下，方程（53）高估了湍动能。基于此，对方程（53）进行重新拟合修正，新的湍动能廓线由下式给出：

图5 湍动能随高度变化规律Fig.5 Variation of turbulent kinetic energy with height

式中，C1和C2为通过拟合获得的常数，其在不同稳定条件下的值由表3 获得。图5（b）给出了修正后的曲线与实测值的对比。从图中可以看出，修正后的曲线与实测值吻合良好。

表3 不同热稳定度下湍动能方程拟合参数Table 3 Fitting parameters for the turbulent kinetic energy equation under different thermal stability degrees

1.4 软件设计及计算方法

1.4.1 软件设计及功能描述

近地面风场特性受地表起伏影响较大。在复杂地形的近地面形成精细湍流结构进而对复杂地形开展CFD 模拟的精度要求很高，导致网格数量庞大、计算耗时长。传统的中尺度大气模拟方法虽然计算速度快，但分辨率在10 km 量级，无法得到精细化的精确风场数据。为了实现精细化大气风场模拟，本文在上述理论基础上，开发了一款具有自主知识产权的风场CFD 模拟软件，能够快速、准确地得到典型地形的风场结果，实现不同大气稳定度条件下的风场分类模拟，可为风资源评估、风电场微观选址和进一步风能预测提供技术支撑，并可节省大量实验成本和建设成本。

自主软件LSWS，即考虑大气稳定度的典型地形风场模拟软件，版本号V1.0。LSWS 软件工作流程及软件构架如图6 所示。软件将大气模式与CFD 结合起来，为典型地形风场模拟提供了大气CFD 模式及模拟方法。软件提供了完备的各种典型地形风场求解方案，入口边界条件可采用多种方案（包括经典的传统大气稳定度廓线、自主修正的Gryning 风廓线等模型），并以实测数据来驱动模式运转，并可以根据使用者的需求来设置入口参数。软件嵌入了适用于我国典型地形、具有不同稳定度修正的廓线模型和相关参数，在CFD 通用方程的基础上增加了热浮力、科氏力等大气运动对流场及湍流的影响模型，具备复杂地形风场大分离、非线性的求解能力，可模拟风场瞬态和稳态流动特征，可有效求解典型地形的包含时间项在内的风场结构。软件包含适用于典型地形风场模拟的网格预处理、几何量处理、求解器与后处理模块，简单易用，且具有跨平台使用能力。

图6 LSWS 软件工作流程及软件架构图Fig.6 Workflow and software architecture of LSWS

1.4.2 计算方法

控制方程的离散采用隐式有限体积法。根据导数项离散采用的不同方法和适用的网格结构类型，离散方法又分为采用局部坐标变换的方法和无坐标变换的方法。本软件的一个特色是给出了无坐标变换方法在斜网格情况下的修正技术，使用该技术分别给出了对流项、黏性项和源项的斜网格修正离散方法。对修正QUICK 类迎风格式，利用边界零体积和同位网格的处理技术，给出了边界无需降阶和外插的处理方案。同时，给出了斜网格修正的边界条件处理方法。软件还开发了每个单元的通量守恒自检功能，可以通过该功能相应地调整网格和计算方法，来保证计算过程中的通量守恒性、收敛性和计算结果的准确性。具体细节可参阅文献[20]。

2 典型地形算例及分析

本文选取锡林浩特和山西南桦山两处典型地形进行大气CFD 模式验证。锡林浩特为平坦的草原地形，山西南桦山为丘陵地形，两者具有代表性，且前期在这两地开展了长期风速测量实验，数据较为完备。图7（a）为锡林浩特平坦草原地形，图7（b）为山西南桦山丘陵地形，图7（c）为南桦山测风点位的山前、山脊和山后分布情况。经过调研和数值实验，本文建立了锡林浩特平坦草原地形计算域模型和山西南桦山丘陵地形计算域模型，如图8 所示。计算域整体采用圆形，便于不同风向的批量计算。在计算域建模方面，将计算域分成平坦区、过渡区和核心区三部分，可使本文提出的风廓线模型既能应用于平坦地形又能应用于复杂地形，同时从计算流体力学角度可使计算域整体满足守恒性，加速收敛，且能避免复杂地形“人造断崖”带来的误差干扰。

图7 锡林浩特和山西南桦山风场位置Fig.7 Locations of the Xilinhot and Nanhuashan wind fields

图8 计算域及网格拓扑Fig.8 Computational domain and grid topology

由于锡林浩特为平坦草原地形，关注区域为6 km范围内的核心区，通过一定的平滑过渡和拉平得到过渡区和平坦区，最终圆形计算域半径为8 km，计算域高度为0.4 km。山西南桦山计算域，考虑到丘陵的相互影响，关注区域在15 km 范围内的核心区，山的当地相对高度为150 m，在山的西侧有1 号和2 号测风点，山的东侧有3 号测风点。同样通过平滑过渡到平坦区，最终山西南桦山圆形计算域半径为18 km，计算域高度为3 km。因两处的实际主导风向均为西南风，因此本文以西南风方向为例进行仿真模拟，使用的入口廓线为1.3.4 节提出的廓线模型。按照计算域中的结构拓扑，利用商业软件进行结构网格划分，首层网格高度为3 m，核心区水平分辨率控制在30 m左右，两款软件采用相同网格和边界条件。

2.1 典型锡林浩特平坦草原地形仿真

分别利用自主软件LSWS 和商业软件Fluent，对锡林浩特典型平坦草原地形开展CFD 仿真。利用锡林浩特实测风数据标定修正廓线模型，选取不稳定、中性和稳定三种稳定度廓线作为入口条件。选取与测风设备点位一致的3 个测点的风速廓线，进行归一化处理，可得到2 个测点位置的风廓线信息，如图9所示。图9（d）为锡林浩特平坦草原地形的海拔高度及测风设备的空间位置。从图9 中可以看出，LSWS和Fluent 计算的结果总体趋势一致，与实测风数据吻合度较高，两软件均能很好地反映出不同稳定度廓线的垂直分布情况，地形对廓线的影响较小。对于不稳定状态（图9a），风速分布斜率较大；中性状态（图9c），斜率接近对数率分布特征；稳定状态（图9e），斜率相对平缓。对于稳定状态风速廓线的计算而言，LSWS相对Fluent，在捕捉廓线廓形上的精度相对更好。

图9 测点位置速度廓线Fig.9 Velocity profiles at the testing locations

测点1 和测点2 在70 m 高度处实测风数据的误差对比如表4 所示。自主软件LSWS 商业软件Fluent误差减小量的计算方式为：Fluent 与实测风速之间的绝对误差E1与LSWS 与实测风速之间的绝对误差E2的差值的百分比，即(E1-E2) × 100%。总体而言，在平坦地形条件下，LSWS 与Fluent 总体精度相当。所有稳定度下，测点1 的误差在1%以内；测点2 在稳定条件下，LSWS 计算精度比Fluent 提升6.33%。

表4 70 m 高度测风数据误差对比Table 4 Error comparison of the wind measurement data at the height of 70 m

2.2 典型山西南桦山丘陵地形仿真

利用自主软件LSWS 和商业软件Fluent，对山西南桦山丘陵地形进行仿真计算对比研究，同时进一步比较LSWS 求解复杂地形风场的能力。同样采用本文的修正廓线模型作为入口边界条件。通过比较CFD 模拟计算结果与实测结果，考察大气CFD 模式及模拟方法对山区复杂地形的影响、风场的空间分布、加减速效应的模拟精度和计算能力。

图10 显示了三种稳定度下、距地面10 m 高度处的速度矢量图。自主软件LSWS 和商业软件Fluent计算的风场分布整体相似。相比之下，Fluent 计算的山区复杂地形迎风侧的风速偏高、背风侧山谷处的风速偏低。

图10 三种稳定度下距地面10 m 高度处速度矢量图Fig.10 Velocity vectors at 10 m above the ground under three stability degrees

自主软件LSWS 和商业软件Fluent 计算得到的三种稳定度下的风场云图如图11 所示。从图中可以看出，对于不稳定状态，两软件的计算结果总体趋势相同；而对于中性状态和稳定状态的风速分布，LSWS 得到的边界层分布更符合真实大气边界层状态，Fluent 计算得到的边界层形成过快、速度梯度较大且对地形的起伏变化不敏感。

图11 三种稳定度下的风场云图Fig.11 Wind field contours under three stability degrees

图12 为3 个典型测点位置的风速廓线分布对比。整体而言，自主软件LSWS 模拟的不同稳定度状态下的风速廓线均与实验值接近，可以较为准确地刻画山地地形动力强迫作用对山前、山脊和山后风速廓线的影响。不同稳定度入口廓线对山地风速廓线的分布影响较小，地形的作用占主导作用。从廓线的总体垂直分布来看，Fluent 在模拟山地风速分布上，速度垂直梯度变化较快，与山地边界层风场特性有一定的偏差。

图12 测点位置风廓线Fig.12 Velocity profiles at the testing locations

1 号、2 号和3 号测风塔在100 m 高度处实测风数据的误差对比如表5 所示。总体而言，在山区复杂地形条件下，自主软件LSWS 的计算精度明显优于商业软件Fluent，所有稳定度下的计算精度均取得提升。对于山前1 号测风塔，在不稳定、中性和稳定条件下，LSWS 与Fluent 计算精度相当。随着气流绕过山脊，对于山脊2 号测风塔，在不稳定、中性和稳定条件下，LSWS 与Fluent 对比，计算精度最大可分别提升5.42%、5.34%、40.77%。就山后3 号测风塔而言，在不稳定、中性和稳定条件下，LSWS 与Fluent 对比，计算精度最大可分别提升38.39%、39.57%和41.80%。由此可见，LSWS 对复杂山区地形影响风场结果的模拟能力有显著的提升。

表5 测风塔100 m 高度测风数据误差对比Table 5 Error comparison of the wind data measured at the height of 100 m

3 结论及建议

本文根据风场实测数据，提出了两种能够反映不同大气稳定度特征的风廓线模型，建立了能够求解不同大气稳定度状态的大气CFD 模式，改进了非中性大气湍流闭合方程及求解方法，提高了可实现性和计算效率，开发了自主可控的CFD 软件，并利用锡林浩特平坦草原地形和山西丘陵地形实测风数据进行了验证。研究得出如下结论：

1）基于实测风数据提出的符合我国地形与风特性的修正Gryning 廓线模型和新的拟合廓线模型，可应用于平坦和复杂地形风场CFD 模拟。

2）通过建立的大气CFD 模式及方法，对锡林浩特平坦草原地形风场进行仿真计算，对比了自主软件LSWS 与商业软件Fluent 的计算精度。在平坦地形下两款软件的计算精度总体相当。对于稳定大气状态，自主软件LSWS 的计算精度较商业软件Fluent 最高可提升6.33%。

3）对山西丘陵地形风场进行仿真计算，总体上自主软件LSWS 的计算精度明显高于商业软件Fluent。在所有大气条件下，自主软件LSWS 较商业软件Fluent 计算精度最高可提升41.80%。特别是对于绕山后流场，本文方法得到的结果更加合理准确，对解决复杂山区风特性评估等问题有积极的意义和工程应用价值。

综合多年CFD 软件开发和商业软件Fluent 使用经验，作者认为本文开发的软件与目前求解域模型和网格更加适配，在计算过程中，还可以查看每一个单元网格的通量守恒性并进行修正，防止了流场的数值扩散，使得计算结果有更好精度。目前自主软件LSWS 还处于起步阶段，操作还不够便捷，需进一步优化改进交互界面，提升用户体验感，进而提升市场竞争力。

致谢：在大气边界层风特性理论、风场实测数据、实测数据处理方法等方面，得到了国家气候中心朱蓉研究员、中国科学院大气物理研究所程雪玲研究员、华北电力大学李莉副教授的支持和指导，在此表示感谢。