一种适用于复杂地形风场流动数值模拟的侧向边界条件设置方法

技术领域

本发明属于计算流体力学技术领域，具体涉及一种适用于复杂地形风场流动数值模拟的侧向边界条件设置方法。

背景技术

我国内陆地形复杂，潜在和现役的大型陆上风电场大多坐落于山地、高原、丘陵、盆地等地带的风资源聚集区域，以计算流体力学(Computational Fluid Dynamics,CFD)为基础的数值模拟方法已成为研究复杂地形风场风资源评估的重要工具。通常，风场计算范围的选取是以地表下垫面为基准面，竖直向上延伸至大气边界层高度(一公里量级)而形成的六面体区域，在数值计算中需要在不同边界面对风速、压强等变量的边界条件进行合理的选定与设置，从而确保流动控制方程在计算区域内的可解性。研究表明，风场计算区域的边界条件设置方法作为数值模拟中的关键因素，直接影响着整个运算过程的收敛性与稳定性，进而决定着数值求解算法的可行性与精确度。

目前，有关复杂地形风场CFD数值模拟边界条件的研究主要集中在入口边界生成与地表边界的处理两个方向上。由于侧向边界条件对数值模拟的空间离散格式、下垫面边缘的不规则程度以及风速-压强耦合算法的敏感型，对其最优化设置一直以来是CFD领域的一个开放性的研究课题。目前，国内外对风场CFD数值模拟中风速与压强条件在侧向边界上的选取与设置，通常是以求解流体力学偏微分(Partial Deferential Equation，PDE)控制方程组(即Navier-Stokes方程组)的两大类边界条件为基础，对第一类Dirichlet边界条件(即待求变量在边界面上的值被直接给定)与第二类Neumann边界条件(即待求变量在边界面外法线的方向导数被指定)进行直接选用、修正或组合使用。

例如，Churchfield M.J.等在侧向边界面采用以Dirichlet边界条件为基础建立的连续性周期(Periodic)边界条件，即两个(一对)侧向边界上所有的变量值对应相等，对山丘与山谷相结合的风场地形上空的大气边界流动进行模拟。然而，周期边界条件要求两个侧向边界中的源边界与目标边界在几何形状与网格划分上必须完全一致，从而实现流动变量从源界面到目标界面的完整映射，这就使得该边界条件仅适用于典型二维地形下流场的模拟验证，无法应用在实际的三维复杂地形中。Peralta C.等在侧向边界面采用以Neumann边界条件为基础的滑移(Slip)边界条件，即流体在侧向边界面上的法向速度为零，切向速度与沿边界面内侧第一层网格的切向速度相等，对坐落于海中的陡崖独立山丘的风场流动进行模拟。然而，滑移边界条件在机理上不允许流体穿过边界面，无法反映现实中流体通量的交换，造成流场模拟在边界区域的不准确性。Balogh M.等在侧向边界面采用Dirichlet与Neumann相结合的对称(Symmetry)边界条件，对缓坡独立山丘周围的风场进行模拟，在一定程度上减少了计算量，但仍旧对实际流动情况做了相应的简化而无法反映气流在边界的真实状态。

综合而言，上述常用的侧向边界条件在计算程序实现上较为简洁，但是在遇到不规则几何边界(如山区、丘陵地带等复杂地形)时，存在着无法适用导致求解发散，或者使用后无法准确模拟边界处风场的实际流动状态，从而对计算域整体的求解准确性造成影响。

发明内容

本发明的目的在于提供了一种适用于复杂地形风场流动数值模拟的侧向边界条件设置方法，该边界条件在求解山区、丘陵等复杂地形地貌上空大气流动时，确保了运算稳定性与收敛性，克服了现有简化型侧向边界条件产生的边界效应对计算精确度的影响。

为了达到上述目的，本发明采用的技术方案是：

本发明提供的一种适用于复杂地形风场流动数值模拟的侧向边界条件设置方法，包括以下步骤：

步骤1，根据目标风电场区域，选定数值模拟的六面体计算域，确定计算域各边界面的位置：将目标区域的复杂地形作为地表边界面，入口和出口边界面垂直于该区域的主风向布置，顶面边界沿水平方向布置并确保最小离地高度为公里量级，侧向边界面平行于主风向布置

步骤2，对步骤1中得到的六面体计算域进行空间网格划分，得到若干个有限体积单元；

步骤3，设定大气流动数值模拟的湍流模式，并根据选定的湍流模式，计算侧向边界面上的各个网格面元对应的有限体积单元体心处，当前求解时间步的前一时间步所对应的流体体积通量

步骤4，根据步骤3中得到的各个网格面元的前一时间步(m-1)所对应的流体体积通量

优选地，步骤3中，所述大气流动数值模拟的湍流模式包括雷诺平均模式和大涡模拟模式。

优选地，若采用雷诺平均模式，通过下式计算各个网格面元的前一时间步所对应的流体体积通量

其中，

优选地，若采用大涡模拟模式，通过下式计算各个网格面元的前一时间步所对应的流体体积通量

其中，下标i为有限体积单元编号；N为目标网格面元周围相邻的网格面元的数目。

优选地，步骤4中，对各个网格面元面心处的风速和压强条件进行设置，具体地：

若

若

式中，

与现有技术方法相比，本发明的有益效果是：

本发明提供的一种适用于复杂地形风场流动数值模拟的侧向边界条件设置方法，能够充分考虑到复杂地形对近地层区域大气流向的局部改变效应，模拟得到的流场特征在侧向边界面附近更加符合实际大气流动状态；

进一步地，本发明提供的速度与压强相耦合的边界条件，能够在同一侧向边界面上模拟出大气入流和出流并存的真实物理现象；

进一步地，本发明提供的速度-压强耦合型侧向边界条件，适用于两种常用的大气湍流模式(雷诺平均和大涡模拟模式)，能够克服崎岖度较高的地形引起的数值求解发散与震荡问题。

附图说明

图1为本发明所涉及的一种适用于复杂地形风场流动数值模拟的侧向边界条件设置方法的流程图；

图2为计算域各边界示意图；

图3为主风向左边侧向边界面中，网格面元及其面心，对应的有限体积单元及其体心示意图；

图4为采用本发明所述侧向边界条件后，主风向右边侧向边界面的横向风速分量模拟云图。

具体实施方式

下面结合附图，对本发明进一步详细说明。

如图1所示，本发明提供的一种适用于复杂地形风场流动数值模拟的侧向边界条件设置方法，包括如下步骤：

步骤1，在目标风场区域确定复杂地形风场流动数值模拟的六面体计算域，并确定各边界面的位置；

根据目标风场水平方向的维度选定数值模拟复杂地形风场的六面体计算域：将目标区域的复杂地形作为地表边界面，入口和出口边界面垂直于该区域的主(盛行)风向布置，顶面边界沿水平方向布置并确保最小离地高度为1公里量级，侧向边界面平行于主(盛行)风向布置，如图2所示。

步骤2，设定大气流动数值模拟的湍流模式，并根据选定的湍流模式，计算步骤1中确定的侧向边界面上的各个网格面元对应的有限体积单元体心处，当前求解时间步(m)的前一时间步(m-1)所对应的流体体积通量

所述大气流动数值模拟的湍流模式根据解析精度通常分为两大类，分别是雷诺平均(Reynolds-Averaged Navier-Stokes，RANS)模式和大涡模拟(Large-eddy Simulation，LES)模式；

计算步骤1中确定的侧向边界面上的各个网格面元对应的有限体积单元体心处，当前求解时间步(m)的前一时间步(m-1)所对应的流体体积通量

(a)若采用RANS模式，则通过下式计算

其中，

(b)若采用LES模式，则通过下式计算

其中，下标i为有限体积单元编号，结合图3所示，i＝0代表目标网格面元所在的有限体积单元，N＝2-4代表目标网格面元周围相邻的网格面元的数目，当目标网格面元位于顶角处时，N＝2；当目标网格面元位于边缘处时，N＝3；当目标网格面元位于其余位置时，N＝4；这里目标网格面元处的

步骤3，根据侧向边界中各面元处求得

(a)若

式中，

(b)若

式中，目标网格面元面心处在当前时间步的切向速度分量

实施例

本发明提供的一种适用于复杂地形风场流动数值模拟的侧向边界条件设置方法，包括以下步骤：

步骤(1)，根据目标风电场区域，选定数值模拟的六面体计算域，确定计算域各边界面的位置；

本实施例以某山区风电场为研究对象，通过Shuttle Radar Topography Mission(SRTM)在线数据库中获得该区域的数字地形高程模型。

数值模拟所需计算域的水平方向维度可根据目标风电场的实际占地面积选取(此算例中取8km×5km)，垂直方向维度选取常用的公里级大气边界层高度(此算例中取1km)。选定的数值模拟计算域中包含一条主山脊(其海拔落差约为400米)以及8座测风塔(以字母“SM”加编号表示)，如图2所示。

计算域的地表边界为目标风场复杂地表下垫面，计算域的入口、出口边界垂直于该区域的主(盛行)风向布置，计算域的侧向边界(两个)平行于主(盛行)风向布置，计算域的顶面边界沿水平方向布置，得到六面体计算域，如图2所示。

步骤(2)，对步骤(1)中得到的六面体计算域进行空间网格划分；

为保证数值模拟的计算精度，采用六面体结构化的网格单元对步骤(2)中得到的六面体计算域进行划分，其网格单元的分辨率原则上须解析出复杂地形在水平和竖直方向的变化特征；

针对本算例地形，所述网格单元的水平方向空间分辨率为15m；垂直方向空间分辨率在近地面进行加密，最底层网格高度为4m，并以tanh的函数关系沿垂向延伸，顶层网格高度为40m。

步骤(3)，根据大气流动数值模拟中采用的湍流模式，分别设置计算域入口、出口、顶部、地表边界的边界条件；

本算例通过选取LES模式对复杂地形风场的大气流动进行数值模拟，基于LES模式的入口边界条件以参考测风塔上最高测量层在主风向内(此算例中为测风塔SM03在57m高处的西南风向)的平均风速为驱动，采用全周期(fully periodic)的前导(Precursor)方法生成；

出口边界条件是采用通用的CFD压力出口边界条件，即压力条件为标准大气压强，风速条件为零梯度变化；

顶面边界条件采用常用的自由滑移边界(Slip)条件；

地表边界条件采用适用于LES模式的

步骤(4)，根据大气流动数值模拟中采用的湍流模式，设置侧向边界面的边界条件：

(4-a)基于大气流动数值模拟中采用的LES模式，对侧向边界面上的各个网格面元对应的有限体积单元体心处，相对于当前求解时间步(m)的前一时间步(m-1)所对应的流体体积通量

其中，

(4-b)通过判定

若

式中，

若

式中，目标网格面元面心处在当前时间步的切向速度分量

步骤(5)，在上述方法中设定的侧向边界条件下，对整个计算域范围大气流动的速度、压强进行CFD数值求解：

本发明中的数值模拟基于开源的CFD场运算和处理软件OpenFOAM(Open SourceField Operation and Manipulation)为平台进行开发，该计算平台以有限体积法(FiniteVolume Method，FVM)为基础，对描述大气湍流运动的不可压缩流体Navier-Stokes方程组进行数值求解，其压强与速度耦合采用PISO(Pressure Implicit with Splitting ofOperators)算法。其中梯度

步骤(6)，数值模拟结果的定性与定量分析：

定性地，在主风向右边侧向边界处的横向风速分量(Uy)模拟云图如图4所示，可以看出，使用本发明所述的速度与压强相耦合的侧向边界条件，能够充分考虑到复杂地形对侧向边界面近地区域大气流向的局部改变效应，同时考虑到由于地球自转产生的地转偏向力对整个风场上空大气流动的宏观转向作用，实现了在侧边界面上同时允许大气入流(Uy为正)和出流(Uy为负)的真实物理现象的模拟。

为进一步验证本发明所述侧向边界条件设置方法的合理性，可将目标区域中8座测风塔位置处的模拟风统计数据(即平均风速与风速标准差)与实测数据进行定量对比，如表1所示。可以看出，模拟得到的风速平均值与波动状况在空间上的变化趋势，与实测数据基本相符；同时，各点位处模拟与实测风数据误差亦在可接受范围(10％)之内。

以上所述，仅为本发明较佳的指导性实施方式，本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明所述的技术范围内，可轻易想到的以该侧向边界条件为核心，对计算平台(如OpenFOAM、Fluent、CFX等)、湍流模式(RANS、LES以及DNS)、以及其它模拟参数的替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权力要求书的保护范围为准。

表1