CN108664705A - 一种基于OpenFOAM的模拟复杂地形地表粗糙度的方法 - Google Patents

Abstract

本发明属于计算风工程领域，并公开了一种基于OpenFOAM的模拟复杂地形地表粗糙度的方法，包含以下步骤：1)建立复杂地形CFD模型，生成精细化网格；2)提取复杂地形CFD模型的地表粗糙度长度数据，将地表粗糙度长度转化为地表覆盖植被高度，为模拟地表覆盖植被对复杂地形风场影响提供准确的数据基础；3)计算地表覆盖植被在不同高度处对风场的阻力系数，然后将地表覆盖植被高度转换为相应的阻力系数；4)使用OpenFOAM添加阻力系数和地表覆盖植被阻力项，从而获得更精确的复杂地形风场特性模拟结果。本发明使用计算机支持的程序语言，考虑了地表粗糙度对复杂地形风场的影响，有利于提升复杂地形风场数值模拟的精确度。

Description

一种基于OpenFOAM的模拟复杂地形地表粗糙度的方法

技术领域

本发明属于计算风工程领域，更具体地，涉及一种模拟复杂地形地表粗糙度的方法。

背景技术

地表粗糙度长度，是指因地表起伏不平或地物本身几何形状的影响，风速廓线上风速为零的位置并不在地表(高度为零处)，而在离地表的一定高度。这一高度也称为空气动力学粗糙度长度，是地表对经过的流体的型态和阻滞力产生影响的一个综合特征参数。地表粗糙度长度是由于地表的植被等产生的，表示了地表的空气动力学特征，反应地表对风速的减弱作用，取决于地面粗糙因素，如粗糙元素的几何形状、大小高度、分布类型和分布密度。研究地表粗糙度长度对分析近地层中气候与环境的变化、风场模拟和风资源评估的准确性具有重要意义。

对于地表粗糙度长度的模拟，一种方法是基于湍流模型的壁面函数法。Blocken等从理论上讨论了传统CFD(Computational Fluid Dynamics,计算流体动力学)设置在模拟粗糙壁面条件风场方面的制约性，并提出了相应的改进方法；Knopp等提出了CFD模拟中适用于模型的粗糙度修正方法。但这些研究没有针对实际地形测试风场分布对地表粗糙度的敏感性，也没有对实际复杂地形中的地表粗糙度长度变化做出适当调整。邓院昌等对整个模拟区域的下垫面进行粗糙度精细化设置来提高模拟精度，但是在数值模拟进行区域分块时均设置了相同的粗糙度，并未最大程度地考虑复杂地形的实际情况。而且使用基于湍流模型的壁面函数法模拟复杂地表粗糙度长度时，无法捕捉粗糙区域内的湍流结构，也没有考虑地表覆盖植被相关参数的影响。

另一种方法是通过在Navier-Stokes动量方程中加入阻力项来模拟，根据Enoki与Ishihara的研究采用如下表达式计算阻力项：

式中，为阻力项，ρ为密度，C_d为阻力系数，为来流风速，为来流风速分量。

采用CFD方法能够模拟大气边界层中的湍流在复杂地形中产生的撞击、分离、环绕等现象；评估复杂地形中的风能分布时具有具有费用低、周期短、便于模拟真实环境、描述风场细节和给出风场定量结果的优点。CFD商业软件简单易用，但是其购买费用昂贵，并且源程序不公开，所提供的二次开发接口有较大局限性，难以满足用户的全部需求。

发明内容

针对现有技术的以上缺陷或改进需求，本发明提供了一种基于OpenFOAM的模拟复杂地形地表粗糙度的方法，该方法使用计算机支持的程序语言，能够基于GIS数据建立复杂地形CFD模型，提取地表粗糙度长度数据，使用OpenFOAM添加阻力系数和地表覆盖植被阻力项，考虑了地表粗糙度对复杂地形风场影响，有利于提升复杂地形风场数值模拟的精确度。

为实现上述目的，按照本发明，提供了一种基于OpenFOAM的模拟复杂地形地表粗糙度的方法，其特征在于，该方法包含以下步骤：

1)在OpenFOAM中基于GIS数据建立复杂地形CFD模型，生成精细化网格；

2)基于GIS数据提取复杂地形CFD模型的地表粗糙度长度数据，然后将一系列分布在等值线上的随机点的地表粗糙度长度转化为地表覆盖植被高度；

3)计算地表覆盖植被在不同高度处对风场的阻力系数，然后将地表覆盖植被高度转换为相应的阻力系数；

4)使用OpenFOAM添加阻力系数和地表覆盖植被阻力项，从而获得复杂地形风场特性模拟结果。

优选地，步骤1)具体包括以下子步骤：

1.1)确定目标的中心经纬度，从GIS地理信息系统中获取目标所在区域的ASTERGDEM30米分辨率大地坐标下的地形数据，使用ArcMap处理GIS数据，从而获得WGS_1984_UTM坐标系下的地形数据网格节点；

1.2)参数化建立三维平坦地形CFD模型并生成精细化网格：复杂地形CFD模型的计算域设置为长方体区域，根据实际的地形地貌找到对计算域有影响的复杂地形，计算域需要包含这些复杂地形，而且采用非结构化三棱柱网格以充分拟合复杂几何边界，其中，在计算域的不同区域采用分辨率不同的网格，复杂地形CFD模型在竖直方向上采用相邻网格尺寸比值为定值的σ网格；

1.3)实现计算域的地形数据网格节点与三维平坦地形CFD模型网格节点之间的数据转换：根据三维平坦地形CFD模型中每一个底面网格数据节点的坐标，找到其在地形数据网格节点中的对应位置，然后使用该位置周围的3个地形数据网格节点的高程，通过三角形共面线性插值得到底面网格节点的高程，同时，在三维平坦地形CFD模型中，竖直方向上其余网格节点的高程会随着底面网格节点高程的变化产生相应变化，最终得到修改高程后的复杂地形CFD模型；

1.4)通过转化器gambitToFoam导入OpenFOAM，将复杂地形CFD模型转换为OpenFOAM可识别的格式。

优选地，步骤2)具体包括以下子步骤：

2.1)从GIS地理信息系统中提取复杂地形CFD模型的地表粗糙度长度数据，地表粗糙度长度数据的分布形式为一系列地表粗糙度长度等值线上的随机点；

2.2)处理地表粗糙度长度数据：判断各随机点是否在复杂地形CFD模型计算域内，如果在计算域内，将地表粗糙度长度转换为地表覆盖植被高度，否则弃之不用；

其中，地表粗糙度长度与地表覆盖植被高度的转换表达式如下：

h＝a·z₀

式中，h为地表覆盖植被高度，a为换算系数，z₀为地表粗糙度长度；

2.3)建立矩形分布的网格节点，使用最近邻点插值法进行赋值，得到矩形分布网格节点的地表覆盖植被高度，其中，最近邻点插值法在搜索时需要设置一个搜索半径，搜索半径不小于矩形分布网格的尺寸，如果搜索半径内没有邻近点，则该节点的地表覆盖植被高度为0。

优选地，步骤3)具体包括以下子步骤：

3.1)根据复杂地形CFD模型底层的每一个三棱柱网格单元的中心坐标(x，y)，找到其在矩形分布的网格中的对应位置，然后使用该中心坐标(x，y)周围的3个地表覆盖植被高度网格节点，通过三角形共面线性插值得到该中心坐标(x，y)处的地表覆盖植被高度；

3.2)从该中心坐标(x，y)处的底层网格单元开始，依次向上，判断每一个三棱柱网格单元距地面的竖直高度z和中心坐标(x，y)处的地表覆盖植被高度之间的关系，如果三棱柱网格单元距地面的竖直高度z不大于中心坐标(x，y)处的地表覆盖植被高度，计算三棱柱网格单元对应的阻力系数，否则阻力系数为0，阻力系数的表达式如下：

C_d＝C_t·a_t

式中，C_d为阻力系数，C_t为抗力系数，a_t为叶面积密度。

优选地，步骤4)具体包括以下子步骤：

4.1)对于复杂地形CFD模型的边界条件，底面采用壁面，入口采用速度入口，出口采用压力出口，其余侧面和顶面采用对称边界条件；

4.2)使用OpenFOAM添加每一个三棱柱网格单元对应的阻力系数；

4.3)在Navier-Stokes动量方程中加入阻力项，考虑了地表粗糙度对风场的影响，使用OpenFOAM计算，从而获得复杂地形风场特性模拟结果，阻力项的表达式如下：

式中，为阻力项，ρ为密度，C_d为阻力系数，为来流风速，为来流风速的速度分量。

总体而言，通过本发明所构思的以上技术方案与现有技术相比，能够取得下列有益效果：

(1)本发明可以实现在OpenFOAM中基于GIS数据建立复杂地形CFD模型，提取地表粗糙度长度数据并转化为地表覆盖植被高度，计算地表覆盖植被在不同高度处对风场的阻力系数，然后将地表覆盖植被高度转换为相应的阻力系数，使用OpenFOAM添加阻力系数和地表覆盖植被阻力项。本发明避开冗杂的手动操作流程，使用计算机支持的程序语言，参数化、程序化和高效化的实现该技术方案。

(2)在复杂地形CFD模拟中，计算域的地表粗糙度长度并非均一不变，相反会随地形地貌呈现出较明显的分布变化。本发明基于GIS数据提取地表粗糙度长度，为模拟地表粗糙度对复杂地形风场的影响提供准的确数据基础。

(3)使用加入阻力项方法模拟地表粗糙度对风场的阻碍作用，不仅考虑了地表覆盖植被中的湍流结构，还考虑了地表覆盖植被相关参数对风场的影响。

附图说明

图1是本发明基于OpenFOAM的模拟复杂地形地表粗糙度的流程图；

图2a和图2b分别是本发明在OpenFOAM中基于GIS数据建立的复杂地形CFD模型和精细化网格分布示意图；

图3a和图3b分别是本发明实施例将地表粗糙度长度数据转换为水平面方向和竖直方向阻力系数的示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

参照各附图，本发明提供的基于OpenFOAM的模拟复杂地形地表粗糙度的方法，其流程图如图1所示。

1、在OpenFOAM中基于GIS数据建立复杂地形CFD建模，生成精细化网格；

1.1)实施例1选取目标的中心经纬度坐标为：经度113.328778°，纬度22.055561°。从GIS地理信息系统中获取目标所在区域的ASTER GDEM30米分辨率大地坐标下的高程数据，使用ArcMap处理GIS数据，从而获得WGS_1984_UTM坐标系下的地形数据网格节点；

1.2)参数化建立三维平坦地形CFD模型并生成精细化网格：三维平坦地形CFD模型的计算域设置为长方体区域，根据实际的地形地貌找到对计算域有影响的复杂地形，计算域需要包含这些复杂地形，而且采用非结构化三棱柱网格以充分拟合复杂几何边界，其中，在计算域的不同区域采用分辨率不同的网格，三维平坦地形CFD模型在竖直方向上采用相邻网格尺寸比值为定值的σ网格；

实施例1的三维平坦地形CFD模型长宽高分别取24km，24km和8km，在关心区域采用细密网格，网格尺寸为20m，远离关心区域采用粗糙网格，竖直方向底面最小网格尺寸为1m，计算域内网格最大生长率为1.2；

1.3)实现计算域的地形数据网格节点与三维平坦地形CFD模型网格节点之间的数据转换：根据三维平坦地形CFD模型中每一个底面网格数据节点的坐标，找到其在地形数据网格节点中的对应位置，然后使用该位置周围的3个地形数据网格节点的高程，通过三角形共面线性插值得到底面网格节点的高程，同时，在三维平坦地形CFD模型中，竖直方向上其余网格节点的高程会随着底面网格节点高程的变化产生相应变化，最终得到修改高程后的复杂地形CFD模型，如图2a、图2b所示；

1.4)通过转化器gambitToFoam导入OpenFOAM，将复杂地形CFD模型转换为OpenFOAM可以识别的格式。

2、基于GIS数据提取复杂地形CFD模型的地表粗糙度长度数据，然后处理地表粗糙度长度数据，将一系列分布在等值线上的随机点的地表粗糙度长度转化为地表覆盖植被高度，为模拟地表粗糙度被对复杂地形风场影响提供准确的数据基础；

2.1)从GIS地理信息系统中提取复杂地形CFD模型的地表粗糙度长度数据，地表粗糙度长度数据的分布形式为一系列地表粗糙度长度等值线上的随机点；

实施例1中，地表粗糙度长度数据格式如下：对于每一条等值线，先有一行信息包括该等值线的地表粗糙度长度和等值线上随机点的个数，再依次罗列这些随机点的平面坐标，地表粗糙度长度取值范围为0～0.5。

2.2)处理地表粗糙度长度数据：判断各随机点是否在复杂地形CFD模型计算域内，如果在计算域内，将地表粗糙度长度转换为地表覆盖植被高度，否则弃之不用；

其中，地表粗糙度长度与地表覆盖植被高度的转换表达式如下：

h＝a·z₀

式中，h为地表覆盖植被高度，a为换算系数，z₀为地表粗糙度长度；

实施例1中将地表粗糙度长度转换为地表覆盖植被高度的换算系数取值7.5。

2.3)建立矩形分布的网格节点，使用最近邻点插值法进行赋值，得到矩形分布网格节点的地表覆盖植被高度，其中，最近邻点插值法在搜索时需要设置一个搜索半径，搜索半径不宜小于矩形分布网格的尺寸，如果搜索半径内没有邻近点，则该节点的地表覆盖植被高度为0；

实施例1中，矩形分布的网格尺寸选取100m，考虑到要覆盖计算域，矩形分布的网格区域长和宽取24.1km，矩形分布的网格节点排列方式为从左至右，从下至上规则分布，设置搜索半径为100m；

3、计算地表覆盖植被在不同高度处对风场的阻力系数，然后将地表覆盖植被高度转换为相应的阻力系数；

3.1)根据复杂地形CFD模型底层的每一个三棱柱网格单元的中心坐标(x，y)，找到其在矩形分布网格中的对应位置，然后使用该中心坐标(x，y)周围的3个地表覆盖植被高度网格节点，通过三角形共面线性插值得到该中心坐标(x，y)处的地表覆盖植被高度；

3.2)从该中心坐标(x，y)处的底层网格单元开始，依次向上，判断每一个三棱柱网格单元距地面的竖直高度z和中心坐标(x，y)处的地表覆盖植被高度之间的关系，如果三棱柱网格单元距地面的竖直高度z不大于中心坐标(x，y)处的地表覆盖植被高度，计算三棱柱网格单元对应的阻力系数，否则阻力系数为0，阻力系数的表达式如下：

C_d＝C_t·a_t

式中，C_d为阻力系数，C_t为抗力系数，a_t为叶面积密度；

实施例1中，抗力系数C_t取值0.4，叶面积密度a_t＝0.04e^{-10(z/h-0.5)(z/h-0.5)}，式中z为三棱柱网格单元距地面的竖直高度，h为地表覆盖植被高度，本发明实施例将地表粗糙度长度数据转换为水平面方向和竖直方向阻力系数的示意图，如图3a、图3b所示。

4、使用OpenFOAM添加阻力系数和地表覆盖植被阻力项，从而获得更精确的复杂地形风场特性模拟结果；

4.1)对于复杂地形CFD模型的边界条件，底面采用壁面，入口采用速度入口，出口采用压力出口，其余侧面和顶面采用对称边界条件；

4.2)使用OpenFOAM添加每一个三棱柱网格单元对应的阻力系数；

4.3)在Navier-Stokes动量方程中加入阻力项，考虑了地表粗糙度对风场的影响，使用OpenFOAM计算，从而获得更精确的复杂地形风场特性模拟结果，阻力项的表达式如下：

式中，为阻力项，ρ为密度，C_d为阻力系数，为来流风速，为来流风速的速度分量。

本发明提供的技术方案基于OpenFOAM的模拟复杂地形地表粗糙度的方法，有利于提升复杂地形风资源评估和风场数值模拟的精确度。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种基于OpenFOAM的模拟复杂地形地表粗糙度的方法，其特征在于，该方法包含以下步骤：

1)在OpenFOAM中基于GIS数据建立复杂地形CFD模型，生成精细化网格；

2)基于GIS数据提取复杂地形CFD模型的地表粗糙度长度数据，然后将一系列分布在等值线上的随机点的地表粗糙度长度转化为地表覆盖植被高度；

3)计算地表覆盖植被在不同高度处对风场的阻力系数，然后将地表覆盖植被高度转换为相应的阻力系数；

4)使用OpenFOAM添加阻力系数和地表覆盖植被阻力项，从而获得复杂地形风场特性模拟结果。

2.根据权利要求1所述的一种基于OpenFOAM的模拟复杂地形地表粗糙度的方法，其特征在于，步骤1)具体包括以下子步骤：

1.1)确定目标的中心经纬度，从GIS地理信息系统中获取目标所在区域的ASTERGDEM30米分辨率大地坐标下的地形数据，使用ArcMap处理GIS数据，从而获得WGS_1984_UTM坐标系下的地形数据网格节点；

1.2)参数化建立三维平坦地形CFD模型并生成精细化网格：复杂地形CFD模型的计算域设置为长方体区域，根据实际的地形地貌找到对计算域有影响的复杂地形，计算域需要包含这些复杂地形，而且采用非结构化三棱柱网格以充分拟合复杂几何边界，其中，在计算域的不同区域采用分辨率不同的网格，复杂地形CFD模型在竖直方向上采用相邻网格尺寸比值为定值的σ网格；

1.3)实现计算域的地形数据网格节点与三维平坦地形CFD模型网格节点之间的数据转换：根据三维平坦地形CFD模型中每一个底面网格数据节点的坐标，找到其在地形数据网格节点中的对应位置，然后使用该位置周围的3个地形数据网格节点的高程，通过三角形共面线性插值得到底面网格节点的高程，同时，在三维平坦地形CFD模型中，竖直方向上其余网格节点的高程会随着底面网格节点高程的变化产生相应变化，最终得到修改高程后的复杂地形CFD模型；

1.4)通过转化器gambitToFoam导入OpenFOAM，将复杂地形CFD模型转换为OpenFOAM可识别的格式。

3.根据权利要求1所述的一种基于OpenFOAM的模拟复杂地形地表粗糙度的方法，其特征在于，步骤2)具体包括以下子步骤：

2.1)从GIS地理信息系统中提取复杂地形CFD模型的地表粗糙度长度数据，地表粗糙度长度数据的分布形式为一系列地表粗糙度长度等值线上的随机点；

2.2)处理地表粗糙度长度数据：判断各随机点是否在复杂地形CFD模型计算域内，如果在计算域内，将地表粗糙度长度转换为地表覆盖植被高度，否则弃之不用；

其中，地表粗糙度长度与地表覆盖植被高度的转换表达式如下：

h＝a·z₀

式中，h为地表覆盖植被高度，a为换算系数，z₀为地表粗糙度长度；

2.3)建立矩形分布的网格节点，使用最近邻点插值法进行赋值，得到矩形分布网格节点的地表覆盖植被高度，其中，最近邻点插值法在搜索时需要设置一个搜索半径，搜索半径不小于矩形分布网格的尺寸，如果搜索半径内没有邻近点，则该节点的地表覆盖植被高度为0。

4.根据权利要求1所述的一种基于OpenFOAM的模拟复杂地形地表粗糙度的方法，其特征在于，步骤3)具体包括以下子步骤：

3.1)根据复杂地形CFD模型底层的每一个三棱柱网格单元的中心坐标(x，y)，找到其在矩形分布的网格中的对应位置，然后使用该中心坐标(x，y)周围的3个地表覆盖植被高度网格节点，通过三角形共面线性插值得到该中心坐标(x，y)处的地表覆盖植被高度；

3.2)从该中心坐标(x，y)处的底层网格单元开始，依次向上，判断每一个三棱柱网格单元距地面的竖直高度z和中心坐标(x，y)处的地表覆盖植被高度之间的关系，如果三棱柱网格单元距地面的竖直高度z不大于中心坐标(x，y)处的地表覆盖植被高度，计算三棱柱网格单元对应的阻力系数，否则阻力系数为0，阻力系数的表达式如下：

C_d＝C_t·a_t

式中，C_d为阻力系数，C_t为抗力系数，a_t为叶面积密度。

5.根据权利要求1所述的一种基于OpenFOAM的模拟复杂地形地表粗糙度的方法，其特征在于，步骤4)具体包括以下子步骤：

4.1)对于复杂地形CFD模型的边界条件，底面采用壁面，入口采用速度入口，出口采用压力出口，其余侧面和顶面采用对称边界条件；

4.2)使用OpenFOAM添加每一个三棱柱网格单元对应的阻力系数；

4.3)在Navier-Stokes动量方程中加入阻力项，考虑了地表粗糙度对风场的影响，使用OpenFOAM计算，从而获得复杂地形风场特性模拟结果，阻力项的表达式如下：

式中，为阻力项，ρ为密度，C_d为阻力系数，为来流风速，为来流风速的速度分量。