CN115908739B

CN115908739B - 一种快速生成复杂地形表面高保真结构网格的方法及***

Info

Publication number: CN115908739B
Application number: CN202310215984.4A
Authority: CN
Inventors: 李泽祥; 冯毅; 李鑫洋; 孙鑫宇; 刘伟毅; 李树民
Original assignee: Chengdu Fluid Power Innovation Center
Current assignee: Chengdu Fluid Power Innovation Center
Priority date: 2023-03-08
Filing date: 2023-03-08
Publication date: 2023-05-05
Anticipated expiration: 2043-03-08
Also published as: CN115908739A

Abstract

本发明涉及一种快速生成复杂地形表面高保真结构网格的方法及***，方法包括：获取具有地形数据的STL格式文件，地形数据包括用于表示地形表面的几何特征的多个三角形面片，且面片通过法线向量以及顶点坐标进行描述，其中，顶点坐标从包括两两垂直的X、Z和Y轴的投影坐标系中获取得到；采用坐标数值提取方法从地形数据中提取出多个坐标数据（X，Y，Z），并将其以离散点的形式进行存储；基于坐标数据生成具有规则排序的地面网格坐标点数组，并从Y方向对地面网格坐标点数组进行拉伸以形成目标网格数组；其中，在XY、ZY平面内可通过点的夹角关系对近地面网格进行正交修正。本发明可以在面对山区丘陵等复杂地形时快速生成高保真结构网格。

Description

一种快速生成复杂地形表面高保真结构网格的方法及***

技术领域

本发明涉及计算流体力学仿真领域，具体涉及一种快速生成复杂地形表面高保真结构网格的方法及***。

背景技术

清洁高效的新能源开发是可持续发展的重要保障。风能作为诸多新能源的代表，具有绿色环保、可再生、储量多、分布广等优点，利用风能发电已经得到飞速发展。在风电项目的开发中，风资源评估是首要的关键一步，直接影响风电项目的效益和风力机的安全。此外，风场计算也是气象预报、飞行安全分析等问题研究的重要基础。

为了能够对某一区域的风资源进行评估或预测，通常会对该区域的地形表面进行仿真模拟，例如，利用网格生成技术建立用于模拟地形表面的地形网格模型。随后，再利用地形网格模型进行相应的仿真或耦合，从而实现对该区域的风场模拟以完成相应的风资源评估、气象预报或者飞行安全分析等等。其中，针对地形表面生成高保真的网格模型是准确地进行风场模拟的关键步骤之一。

在现有的网格生成技术中，一般有生成结构网格、非结构网格两种方式。结构网格在实际应用上往往只适用于形状规则的图形，而地形表面（尤其是山地、丘陵地形等）大都存在不同程度的起伏变化，即地形表面属于复杂、不规整的几何外形，因此，针对地形表面生成贴体的结构网格是很困难的。因而在现有的风场模拟等应用中通常选用的是生成非结构网格的方式。例如，申请公布号为CN108563867A的中国发明专利申请，其公开了一种基于OpenFOAM实现WRF和CFD耦合模拟风场的方法，该方法中地形CFD模型采用的则是非结构的三棱柱网格。其中，非结构网格一方面在计算时需要较大的内存，这对计算机的算力提出了更高的要求，另一方面由于非结构网格对边界的恢复比较困难，很可能造成网格生成的失败。

因此，亟需一种可以快速生成高保真地面网格的方法。

发明内容

本发明的目的在于提供一种快速生成复杂地形表面高保真结构网格的方法及***，部分地解决或缓解现有技术中的上述不足，能够快速、准确地生成贴体的地形表面的结构网格。

为了解决上述所提到的技术问题，本发明具体采用以下技术方案：

提供一种快速生成复杂地形表面高保真结构网格的方法，包括步骤：

S101获取具有地形数据的STL格式文件，所述地形数据包括：用于表示所述复杂地形表面的几何特征的多个三角形面片，所述三角形面片通过法线向量以及顶点坐标进行描述，其中，所述顶点坐标从投影坐标系中获取得到，所述投影坐标系包括：X方向和Z方向，以及与所述X方向、Z方向相互垂直的垂向Y方向；

S102采用坐标数值提取方法从所述地形数据中的目标范围内提取出多个坐标数据，并将所述多个坐标数据以离散点的形式进行存储，其中，所述坐标数据表示为（X，Y，Z）；

S103基于所述坐标数据生成具有规则排序的地面第一层网格数组，并从所述Y方向对所述地面第一层网格数组进行拉伸以创建形成目标网格数组；

其中，S103包括步骤：

S31确定目标计算域，并获取所述目标计算域内的所有离散点的坐标数据，利用数组创建方法根据所获取的坐标数据创建得到一维数组Xt、Yt、Zt，其中，所述一维数组Xt、Yt、Zt分别为在X、Y、Z方向上的坐标数值的集合，每个方向上有NX×NZ个坐标数据；

S32按照一维数组Xt中的坐标数值大小对所述一维数组Xt进行第一次排序，并同步变更一维数组Yt和Zt中对应坐标数值Y、Z的顺序，以对应地获得一维数组Xn、Yn、Zn；

S33采用数组变换方法将完成第一次排序的所述一维数组Xn、Yn、Zn变换为尺寸为[NX][NY][NZ]的三维数组Xc、Yc、Zc，其中，所述三维数组Xc、Yc、Zc包括在Y方向上的地面的第一层网格数据；

S34按照三维数组Zc中的坐标数值沿着Z方向的大小对所述三维数组Zc进行第二次排序，并同步变更三维数组Xc和Yc中对应坐标数值的位置，使得所述三维数组Xc、Yc、Zc形成在Y方向上的地面第一层网格数组；

S35根据预设参数计算所述目标计算域内从地面至所述目标计算域的顶部沿Y方向的网格分布系数，并通过所述网格分布系数对所述地面第一层网格数组进行拉伸，从而拉伸创建得到在所述地面第一层网格之上且尺寸为[NX][1][NZ]的NY-1层网格数组，所述地面第一层网格数组和所述NY-1层网格数组共同形成计算域网格数组；

S36按照目标分辨率对所述计算域网格数组进行采样以获取到目标网格点的坐标数值X、Z；

S37基于所述目标网格点的所述坐标数值输出所述目标网格数组。

在一些实施例中，在S35之前，还包括步骤：

在层网格中，从第二层网格开始，通过下层网格与相邻上层网格之间的夹角关系，分别在至少一个平面内对至少一个所述上层网格中的至少一个离散点的坐标进行修正；

其中，在XY平面内，对网格编号为的离散点D的坐标数值X进行修正的修正量为：

；

所述离散点D修正后的坐标数值X'为：

；

其中，和分别为XY平面内的两个特征角∠CBD、∠ABD，离散点A₁的坐标为：,网格编号为，离散点B₁的坐标为：，网格编号为，离散点C₁的坐标为：，网格编号为，离散点D的坐标为：，网格编号为；

在ZY平面内，对网格编号为的离散点D的坐标数值Z的修正量为：；

所述离散点D修正后的坐标数值Z'为：

；

其中，和分别为ZY平面内的两个特征角∠CBD、∠ABD，离散点A₂的坐标为：，网格编号为，离散点B₂的坐标为：，网格编号为，离散点C₂的坐标为：，网格编号为，离散点D的坐标为：，网格编号为。

在一些实施例中，在至少一个平面内对至少一个所述上层网格中的至少一个离散点的坐标进行修正之前，还包括步骤：

计算所述层网格中每一个已修正的离散点的X坐标值在Y方向上的修正量的总和；

当所述总和小于预设的第一阈值时，则继续在当前的所述上层网格中进行修正，否则，停止在所述上层网格中进行修正。

在一些实施例中，计算所述层网格中每一个已修正离散点的Z坐标值在Y方向上的修正量的总和；

当所述总和小于预设的第二阈值时，则继续在当前的所述上层网格中进行修正，否则，停止在所述上层网格中进行修正。

在一些实施例中，在至少一个平面内对至少一个所述上层网格中的至少一个离散点的坐标进行修正的步骤之前，还包括步骤：

计算当前上层网格中每一个已修正的所述离散点的X坐标值在Y方向上的修正值总和；

当所述总和小于预设的第一层修正阈值时，则继续在当前的所述上层网格中进行修正，否则，停止在所述上层网格中进行修正。

在一些实施例中，计算当前上层网格中每一个已修正的所述离散点的Z坐标值在Y方向上的修正值总和；

当所述总和小于预设的第二层修正阈值时，则继续在当前的所述上层网格中进行修正，否则，停止在所述上层网格中进行修正。

在一些实施例中，还包括步骤：

S104设定边界条件，将所述目标网格数组导出为用于模拟所述复杂地形表面的目标结构网格。

在一些实施例中，S31包括：

读取用户预设的顶点坐标数据，并基于所述顶点坐标数据确定所述目标计算域，其中，顶点包括：第一顶点（Xmin、Zmin）、第二顶点（Xmax、Zmax）；

基于所述顶点坐标数据判断所述离散点是否位于所述目标计算域内，并将位于所述目标计算域内部的所述离散点识别为有效坐标点；

利用数组创建方法根据所述有效坐标点的坐标数据创建得到一维数组Xt、Yt、Zt。

在一些实施例中，NX为一维数组Xn中具有同一数值大小的坐标数值个数，NZ为沿Z方向的离散点的个数；NZ为沿Z方向的离散点的个数；其中，NZ=N/NX，N为所述有效坐标点的个数。

在一些实施例中，S35包括：

读取用户预设的垂向网格点坐标增长率调节因子Fac、目标计算域的顶部高度Ytop；

基于Fac、Ytop使用双曲正切函数tanh构造元素范围处于(0,1]区间内的数组H(j)，并计算相邻垂向离散点之间的间距系数dH(j)；

通过逐点累加间距的方式计算得到所述目标计算域内所有离散点的垂向坐标Yc。

在一些实施例中，S36中采样的方法包括：邻点插值法，和/或间隔取点法。

本发明还对应于上述方法提供了一种快速生成复杂地形表面高保真结构网格的***，包括：

数据获取模块，被配置为用于获取具有地形数据的STL格式文件，所述地形数据包括：用于表示所述复杂地形表面的几何特征的多个三角形面片，所述三角形面片通过法线向量以及顶点坐标进行描述，其中，所述顶点坐标从投影坐标系中获取得到，所述投影坐标系包括：X方向和Z方向，以及与所述X方向、Z方向相互垂直的垂向Y方向；

数据提取模块，被配置为用于采用坐标数值提取方法从所述地形数据中的目标范围内提取出多个坐标数据，并将所述多个坐标数据以离散点的形式进行存储，其中，所述坐标数据表示为（X，Y，Z）；

网格生成模块，被配置为用于基于所述坐标数据生成具有规则排序的地面第一层网格数组，并从所述Y方向对所述地面第一层网格数组进行拉伸以创建形成目标网格数组；

其中，网格生成模块包括：

计算域确定单元，被配置为用于确定目标计算域，并获取所述目标计算域内的所有离散点的坐标数据，利用数组创建方法根据所获取的坐标数据创建得到一维数组Xt、Yt、Zt，其中，所述一维数组Xt、Yt、Zt分别为在X、Y、Z方向上的坐标数值的集合，每个方向上有NX×NZ个坐标数据；

一次排序单元，被配置为用于按照一维数组Xt中的坐标数值大小对所述一维数组Xt进行第一次排序，并同步变更一维数组Yt和Zt中对应坐标数值Y、Z的顺序，以对应地获得一维数组Xn、Yn、Zn；

数组变换单元，被配置为用于采用数组变换方法将完成第一次排序的所述一维数组Xn、Yn、Zn变换为尺寸为[NX][NY][NZ]的三维数组Xc、Yc、Zc，其中，所述三维数组Xc、Yc、Zc包括在Y方向上的地面的第一层网格数据；

二次排序单元，被配置为用于按照三维数组Zc中的坐标数值沿着Z方向的大小对所述三维数组Zc进行第二次排序，并同步变更三维数组Xc和Yc中对应坐标数值的位置，使得所述三维数组Xc、Yc、Zc形成在Y方向上的地面第一层网格数组；

拉伸单元，被配置为用于根据预设参数计算所述目标计算域内从地面至所述目标计算域的顶部沿Y方向的网格分布系数，并通过所述网格分布系数对所述地面第一层网格数组进行拉伸，从而拉伸创建得到在所述地面第一层网格之上且尺寸为[NX][1][NZ]的NY-1层网格数组，所述地面第一层网格数组和所述NY-1层网格数组共同形成计算域网格数组；

采样单元，被配置为用于按照目标分辨率对所述计算域网格数组进行采样以获取到目标网格点的坐标数值X、Z；

输出单元，被配置为用于基于所述目标网格点的所述坐标数值输出所述目标网格数组。

有益技术效果：

本发明首先从原始地形数据（例如，从GIS***中导出的矢量地图数据）入手，在综合考量丘陵、山地等复杂地形表面的几何特性（即地势起伏大，整个地形表面在划分结构网格过程中可被视为处于正六面体的下表面）以及所获取的“投影坐标系”的坐标数据格式特点后，引入二次排序、近地面拉伸等关键手段用于生成近地面关键区域的结构网格。其中，正由于本发明提出了二次排序并配合近地面拉伸的网格生成方式，所生成的结构网格更便于在后续进行精度上的调整（也即更利于对近期面区域的网格进行正交修正）。

进一步地，针对二次排序、近地面拉伸处理之后所生成的一种新的结构网格，本发明还提出了利用相邻微面元的夹角关系对离散点（也被称为“网格点”）的水平坐标进行修正的方式，将对三维空间网格结构的修正简化在两个平面（即XY、ZY平面）上进行，大大地减小了修正过程的数据计算量，进而显著地降低了网格修正过程（尤其是针对大范围的复杂地形网格进行精度修正时）对计算机算力的需求。

因此，本发明所提供的网格生成方法在二次排序、近地面拉伸以及正交修正方式等多手段的协同作用下，能够快速地针对近地面关键区域生成高保真的结构网格，从而能够快速、便捷地对近地面关键区域进行准确的风资源评估。这种准确的风资源评估不仅对复杂山区、丘陵地形的风力发电设备或高层建筑（尤其是特大型桥梁等）的安全作业、建设有着关键性的影响，对于运输机或无人机实施空投物资也有重要意义。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作一简单地介绍。在所有附图中，类似的元件或部分一般由类似的附图标记标识。附图中，各元件或部分并不一定按照实际的比例绘制。显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其它的附图。

图1a为本发明一示例性实施例中的快速生成结构网格的流程示意图；

图1b为本发明一示例性实施例中通过二次排序、近地面拉伸生成结构网格的流程示意图；

图2为本发明一具体实施例中的快速生成结构网格的流程示意图；

图3为本发明一具体实施例中网格划分过程的流程示意图；

图4为复杂丘陵地形示意图；

图5本发明一示例性实施例中进行正交修正的特征角度(α-β)/2的示意图；

图6为本发明一具体实施例中正交修正前网格的最大偏斜角度分布的示意图；

图7为本发明一具体实施例中正交修正后网格的最大偏斜角度分布的示意图；

图8为本发明一具体实施例中网格划分过程生成的复杂地形表面网格及垂向网格分布的示意图；

图9为本发明一具体实施例中网格划分过程生成的复杂地形网格计算域整体视图的示意图；

图10为本发明一示例性实施例中的快速生成复杂地形表面高保真结构网格的***的模块结构示意图。

其中，10为数据获取模块，20为数据提取模块，30为网格生成模块。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本文中，使用用于表示元件的诸如“模块”、“部件”或“单元”的后缀仅为了有利于本发明的说明，其本身没有特定的意义。因此，“模块”、“部件”或“单元”可以混合地使用。

本文中，术语“上”、“下”、“内”、“外”“前”、“后”、“一端”、“另一端”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。此外，术语“第一”、“第二”仅用于描述目的，而不能理解为指示或暗示相对重要性。

本文中，除非另有明确的规定和限定，术语“安装”、“设置有”、“连接”等，应做广义理解，例如“连接”，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以具体情况理解上述术语在本发明中的具体含义。

本文中“和／或”包括任何和所有一个或多个列出的相关项的组合。

本文中“多个”意指两个或两个以上，即其包含两个、三个、四个、五个等。

需要说明的是，在本文中，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者装置不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者装置所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括该要素的过程、方法、物品或者装置中还存在另外的相同要素。

如在本说明书中使用的，术语“大约”，典型地表示为所述值的+/-5%，更典型的是所述值的+/-4%，更典型的是所述值的+/-3 %，更典型的是所述值的+/-2 %，甚至更典型的是所述值的+/-1 %，甚至更典型的是所述值的+/-0.5%。

在本说明书中，某些实施方式可能以一种处于某个范围的格式公开。应该理解，这种“处于某个范围”的描述仅仅是为了方便和简洁，且不应该被解释为对所公开范围的僵化限制。因此，范围的描述应该被认为是已经具体地公开了所有可能的子范围以及在此范围内的独立数字值。例如，范围1〜6的描述应该被看作已经具体地公开了子范围如从1到3，从1到4，从1到5，从2到4，从2到6，从3到6等，以及此范围内的单独数字，例如1，2，3，4，5和6。无论该范围的广度如何，均适用以上规则。

本文中，“近地面层（near surface layer）”又称地面边界层（surface boundarylayer），是大气边界层中稳定存在的，最接近下地面的部分。其中，近地面层的高度通常在50至100 m左右，可以通过观测的风速和气温廓线确定，也可通过相似性理论得到参数化的高度，例如Monin-Obukhov长度尺度。

本文中，“近地面层网格”也称为“近地面关键区域网格”，其为用于描述或模拟地形表面至近地面层之间的空间区域的多层结构网格。或者，“近地面层网格”也可为用于描述由用户（如工程人员）自定义的近地面区域的多层结构网格。其中，近地面区域即紧邻地形表面的空间区域，其可以由用户根据实际工程应用需求或工作经验进行适应性的设定。

计算网格按网格点之间的连接关系可分为结构网格（structured grid）和非结构网格（unstructured grid）。其中，结构网格很容易地实现区域的边界拟合，适用于流体和表面应力集中等方面的计算，因此结构网格在CFD计算过程中相对于非结构网格具有显著的优势。但是，结构网格的网格点之间的邻接是有序、规则的，单元是二维的四边形、三维的六面体。结构网格的结构性、有序性会在一定程度上限制其对复杂几何构型/几何外形的适应能力，因此在面临复杂几何外形（如地形表面）的网格问题时其网格生成较困难，网格生成的人工工作量比非结构网格显著增多，对计算机算力要求也相对较高。目前，对于常见领域的复杂几何外形的流动问题，本领域技术人员也有尝试采用分区网格、重叠网格技术等方法克服其几何适应能力差的缺陷，但这又会带来分区网格各区间隐式边界条件难以实现、重叠网格各区间寻找变量传递对应点及变量插值不守恒性等一系列问题。因此，如何能够快速地生成高保真的结构网格是本发明所要解决的关键技术问题之一。

实施例一

本发明第一方面提供了一种快速生成复杂地形表面的高保真结构网格的方法，如图1a和图1b所示，所述方法包括步骤：

步骤S101：获取具有地形数据的STL格式文件，所述地形数据包括：用于表示所述复杂地形表面的几何特征的多个三角形面片，所述三角形面片通过法线向量以及顶点坐标进行描述，其中，所述顶点坐标从投影坐标系中获取得到，所述投影坐标系包括：X方向和Z方向，以及与所述X方向、Z方向相互垂直的垂向Y方向；

在一些具体实施例中，可通过GIS***（geographical information system, 地理信息***）获取复杂地形表面的矢量地图数据，并将其导出为STL格式。其中，STL格式文件可将复杂的实体地形表面（如图4所示的复杂丘陵地形）通过大量的三角形面片进行表示，每个面片由一个垂直方向（即法线向量）和代表三角形顶点（角）的三个点描述。本实施例优选格式简单、只描述三维地形的几何信息的STL格式文件以实现后续准确、快速地结构网格的划分。

步骤S102：采用坐标数值提取方法从所述地形数据中提取出多个坐标数据，并将所述多个坐标数据以离散点的形式进行存储，其中，所述坐标数据表示为（X，Y，Z）。

在一些具体实施例中，将STL格式文件导出为DAT格式的离散点（也即地形坐标点）。本实施例中使用Tecplot软件将STL格式文件导出为DAT格式的离散点数据，借助Tecplot软件，可以便捷地将刻画地形表面的离散点坐标保存为DAT格式数据文件，从而避免出现重复坐标数据（如相邻两个三角形共用的顶点坐标被重复读取）。

步骤S103：基于所述坐标数据生成具有规则排序的地面第一层网格数组，并从所述Y方向对所述地面第一层网格数组进行拉伸以创建形成目标网格数组。

在一些具体实施例中，优选地根据二次排序、近地面坐标拉伸生成近地面区域的目标网格数组，并优选地将其输出为PLOT3D格式。

优选地，本实施例中S103包括步骤：

S33采用数组变换方法将完成第一次排序的所述一维数组Xn、Yn、Zn变换为尺寸为[NX][NY][NZ]的三维数组Xc、Yc、Zc，其中，所述三维数组Xc、Yc、Zc包括在Y方向上的地面的第一层网格数据（也即近地面的第一层网格数据）；

其中，所生成的NY-1层网格中的每一层网格的维度/尺寸为[NX][1][NZ]，网格点总数是NX*1*NZ，整个目标计算域中的网格点总数为NX*NY*NZ。

本实施例中将原始坐标数值转换为了三维数组的格式，从而使得网格点（或者说离散点）更适用于表达原始地形的三维空间区域内的地形坐标点分布特征，进而有助于后续实现平面上的正交修正。

进一步地，在一些实施例中，为了保证近地面层网格的高保真性，优选地对近地面层网格进行正交修正，以提高近地面层网格的质量，在S36之前，还包括步骤：

在所述层网格（优选为近地面层网格）中，从第二层网格开始，通过下层网格与相邻上层网格（也即当前待修正的层网格）之间的夹角关系，分别在至少一个平面内对至少一个所述上层网格中的至少一个离散点的坐标进行修正。

优选地，在一些实施例中，采用正交修正的方式对离散点的水平坐标进行修正，在XY平面内，对网格编号为的离散点D的坐标数值X进行修正的修正量为：

；

所述离散点D修正后的坐标数值X'为：

；

其中，和分别为XY平面内的两个特征角∠CBD、∠ABD，离散点A₁的坐标为：,网格编号为，离散点B₁的坐标为：，网格编号为，离散点C₁的坐标为：，网格编号为，离散点D的坐标为：，网格编号为。

可以理解的是，在ZY平面可以采用相同的正交修正对离散点的另一个水平坐标进行修正。

例如，在ZY平面内，对网格编号为的离散点D的坐标数值Z的修正量为：；

所述离散点D修正后的坐标数值Z'为：

；

可以理解的是，本发明实施例中的网格编号是基于二次排序的结果而确定的。其中，网格编号反应了离散点的坐标数值在二次排序之后的存储位置分布。例如，i、j、k分别表示坐标数值在X、Y、Z方向上的排序编号（或者说，存储位置分布）。

进一步地，在一些实施例中，在逐层向上对水平坐标进行修正的过程中，设置调整的修正阈值范围并使修正量逐渐减小，以避免在对近地面层网格进行修正的过程中，出现负体积或网格尺寸变化过快。

本发明实施例中所关注的风资源分析主要应用于风力发电等工程领域，风力发电机的轮毂高度在100m量级左右，而大气边界层近地面风速分布的准确模拟对评估风力发电等工程应用很关键。除了风力发电，城市建筑、桥梁等的高度也大都在几十米至几百米高度范围内，对高层建筑、桥梁进行空气动力学分析也需要保证近地面关键网格的保真性。本实施例通过对关键近地面区域进行网格修正，其中，在逐层向上修正的过程中，设置调整的阈值（如第一阈值等）并使修正量逐渐减小，既能够保证靠近地面区域的网格有高保真性，同时能够避免出现负体积或网格尺寸变化过快等问题。

例如，在一些实施例中，在至少一个平面内对至少一个所述上层网格中的至少一个离散点的坐标进行修正之前，还包括步骤：

计算所述近地面层网格中每一个修正后的离散点的X坐标值在Y方向上的已经修正后的修正量的总和；

当所述总和小于预设的第一阈值（可由用户根据工程应用需求自行设定）时，则继续在当前的所述上层网格中进行修正，否则，停止在所述上层网格中进行修正。

具体地，在一些实施例中，计算每一个水平坐标点位置处在垂直方向上的修正量的总和，即计算网格点修正量的总和，代表或，n是当前层网格在垂直方向（Y方向）的排列编号。

在XY平面内进行修正时，计算在所述近地面层网格中每一个水平坐标点（固定i、k指标）位置处垂直方向上已经修正过的所述修正量的总和。当所述总和小于预设的第一阈值，即时，则继续在所述当前平面内对所述坐标进行修正，修正量为，否则，停止在所述当前平面内对所述坐标进行修正。

在ZY平面内进行修正时，计算在所述近地面层网格中每一个水平坐标点（固定i、k指标）位置处垂直方向上已经修正过的所述修正量的总和。当所述总和小于预设的第二阈值，即时，则继续在所述当前平面内对所述坐标进行修正，修正量为，否则，停止在所述当前平面内对所述坐标进行修正。

又例如，在一些实施例中，在至少一个平面内对至少一个所述上层网格中的至少一个离散点的坐标进行修正之前，还包括步骤：

计算所述近地面层网格中每一个已修正的离散点的Z坐标值在Y方向上的已经修正后的修正量的总和；当所述总和小于预设的第二阈值时，则继续在当前的所述上层网格中进行修正，否则，停止在所述上层网格中进行修正。

或者，在另一些实施例中，可以对每一层的层网格设定相应的修正阈值，从而对每一层网格的修正量进行针对性的调整。具体地，本实施例中，在至少一个平面内对至少一个所述上层网格中的至少一个离散点的坐标进行修正的步骤之前，还包括步骤：

计算当前上层网格中每一个已修正的所述离散点的X坐标值在Y方向上的修正值总和；当所述总和小于预设的第一层修正阈值时，则继续在当前的所述上层网格中进行修正，否则，停止在所述上层网格中进行修正。同样地，在对Z坐标值进行修正时，还包括：计算当前上层网格中每一个所述离散点的Z坐标值在Y方向上的修正值总和；当所述总和小于预设的第二层修正阈值时，则继续在当前的所述上层网格中进行修正，否则，停止在所述上层网格中进行修正。

优选地，在一些实施例中，层网格的修正阈值从下到上逐渐减小，从而避免出现负体积或网格尺寸变化过快等问题。

进一步地，在一些实施例中，该方法还包括步骤：设定边界条件，将所述目标网格数组导出为用于模拟所述地形表面的目标结构网格。

进一步地，在一些实施例中，S31包括：

读取用户预设的顶点坐标数据，并基于所述顶点坐标数据确定所述目标计算域，其中，顶点包括：第一顶点（Xmin、Zmin）、第二顶点（Xmax、Zmax）；基于所述顶点坐标数据判断所述离散点是否位于所述目标计算域内，并将位于所述目标计算域内部的所述离散点识别为有效坐标点；

其中，NX为一维数组Xn中具有同一数值大小的坐标数值个数，NZ=N/NX，N为所述有效坐标点的个数。

进一步地，在一些实施例中，S35包括：

读取用户预设的垂向网格点坐标增长率调节因子Fac、目标计算域的顶部Ytop（具体地，目标计算域为地面以上的空气流动区域，如目标计算域为近地面的空气流动区域）；

基于Fac、Ytop使用双曲正切函数tanh构造元素范围处于(0,1]区间内的数组H(j)，并计算相邻垂向离散点之间的间距系数dH(j)（也即网格分布系数）；

参见图2-图9所示，下文将结合具体实施例对步骤S103进行示例性说明。另需要说明的是下文中给出的步骤S103实现方式仅作为示例性的说明，并不代表本申请实施例的全部实现方式。

优选地，本实施例中利用Fortran语言程序实现自动化网格生成，从而避免或减少工程人员手动划分网格的工作。其中，Fortran语言程序示例的示例如下：

S31包括步骤：读取DAT格式的离散点坐标数据，同时读取工程人员指定的计算域边界顶点的坐标值Xmin、Xmax、Zmin、Zmax，以及目标计算域的顶部高度Ytop、垂向网格点坐标增长率调节因子Fac等预设参数。其中，Xmin和Xmax代表了水平方向X的最小和最大坐标，Zmin和Zmax代表了另一个水平方向Z的最小和最大坐标，地形网格将在(Xmin,Zmin)和(Xmax,Zmax)构成的矩形区域（即目标计算域）内生成；

判断离散点是否位于目标计算域边界范围内，将目标计算域边界范围内的点计入一维数组Xt、Yt、Zt中，共计N个有效坐标点；其中，逐点判断DAT格式数据文件中的离散点是否位于(Xmin,Zmin)和(Xmax,Zmax)构成的矩形区域内，Fortran语言程序示例为：

N=0

open(10,file='terrain.dat')

do while(not(eof(10)))

read(10,*) Xtemp,Ytemp,Ztemp

if(Xmax>Xtemp .and. Xmin<Xtemp .and. Zmax>Ztemp .and. Zmin<Ztemp)then

Xt(N)=Xtemp

Yt(N)=Ytemp

Zt(N)=Ztemp

N=N+1

endif

enddo

S32包括步骤：按照Xt数组元素（即坐标数值X）的大小进行排序，并同步变更Yt和Zt数组中对应元素（即坐标数值Y、Z）的顺序，获得数组Xn、Yn、Zn；其中，按照Xt数组元素的大小对Xt、Yt、Zt进行同步排序，Fortran语言程序示例为：

do i=1,N-1

do m=i+1,N

if(Xt(i)>Xt(m)) then

Xtemp=Xt(i)

Xt(i)=Xt(m)

Xt(m)=Xtemp

Ytemp=Yt(i)

Yt(i)= Yt(m)

Yt(m)=Ytemp

Ztemp=Zt(i)

Zt(i)=Zt(m)

Zt(m)=Ztemp

endif

enddo

Xn=Xt

Yn=Yt

Zn=Zt

其中，根据STL格式地形数据的特点，Xn数组中具有同一坐标数值的元素个数即为目标计算域边界范围内沿水平方向X的坐标点个数NX，沿水平方向Z的坐标点个数为NZ=N/NX；

判断在(Xmin,Zmin)和(Xmax,Zmax)构成的矩形区域内沿着X和Z方向的网格点数，Fortran语言程序示例为：

NZ=1

do m=2,N

if(Xn(m)==Xn(1)) then

NZ=NZ+1

endif

enddo

NX=N/NZ

S33、S34包括步骤：将已经按照X方向数值大小同步排序的一维数组Xn、Yn、Zn按照对应元素位置同步变更的方式变形为尺寸[NX][NY][NZ]的三维数组Xc、Yc、Zc，并按照三维数组Zc中的坐标数值大小对三维数组Xc和Yc进行第二次排序，本实施例中第二次排序所采用的Fortran语言程序示例为：

do i=0,NX-1

Xc(i,0,0:NZ-1)=Xn(i*NZ+1:i*NZ+NZ)

Yc(i,0,0:NZ-1)=Yn(i*NZ+1:i*NZ+NZ)

Zc(i,0,0:NZ-1)=Zn(i*NZ+1:i*NZ+NZ)

enddo

do i=0,NX-1

do k=0,NZ-1-1

do m=k+1,NZ-1

if(Zc(i,0,k).gt.Zc(i,0,m)) then

Ztemp=Zc(i,0,k)

Zc(i,0,k)=Zc(i,0,m)

Zc(i,0,m)=Ztemp

Ytemp=Yc(i,0,k)

Yc(i,0,k)=Yc(i,0,m)

Yc(i,0,m)=Ytemp

Xtemp=Xc(i,0,k)

Xc(i,0,k)=Xc(i,0,m)

Xc(i,0,m)=Xtemp

endif

enddo

例如，一维数组到三维数组的变化示例为：一维数组Xn(1:N)，元素个数是N，变换后是三维数组Xc(0:NX-1, 0, 0:NZ-1)，元素个数是NX*1*NZ （元素从0开始编号），N=NX*1*NZ。Xn和Xc存储的都是所有地面坐标点的X坐标。

S35包括步骤：读取垂向网格点坐标增长率调节因子Fac和目标计算域的顶部高度Ytop（Fac、Ytop可以由工程人员基于地形表面的起伏状况以及对网格的精度要求进行适应性地设定），使用双曲正切函数tanh构造元素范围处于(0,1]区间内的数组H(j)，并计算相邻垂向坐标点之间的间距系数dH(j)，通过逐点累加间距的方式计算得到计算域内所有网格点的垂向坐标Yc；

获得计算域内所有网格点的垂向坐标数组Yc(0:NX-1,0:NY-1,0:NZ-1)，Fortran语言程序示例为：

do j=0,NY-1

H(j)=tanh((j-NY+1)/(NY*1.0)*Fac)+1.0

enddo

do j=1,NY-1

dH(j)=H(j)-H(j-1)

enddo

do k=0,NZ-1

do i=0,NX-1

d(i,k)= Ytop-Yc(i,0,k)

enddo

do i=0,NX-1

do k=0,NZ-1

do j=1,NY-1

Yc(i,j,k)=Yc(i,j-1,k)+dH(j)*d(i,k)

enddo

此时，由于网格是从第一层网格（地面网格）开始垂直向上拉伸，垂向每一层网格的水平方向坐标数组Xc(0:NX-1,j,0:NZ-1)和Zc(0:NX-1,j,0:NZ-1)与第一层网格的水平方向坐标数值一致，即为Xc(0:NX-1,0,0:NZ-1)和Zc(0:NX-1,0,0:NZ-1)。

为了进一步地保证网格的近地面关键区域具有高保真性，在S34之后，还将对网格点的坐标进行修正，从第二层网格开始，依据相邻低层网格中相邻微面元的边与边之间的夹角关系（如图5所示），分别计算XY和ZY平面内相邻下层网格的相邻边的角平分线与Y轴的夹角，根据夹角值分别对该层网格的网格点的水平坐标X和Z坐标进行调整，得到修正后的计算域内所有网格点的水平方向坐标数组Xc(0:NX-1,0:NY-1,0:NZ-1)和Zc(0:NX-1,0:NY-1,0:NZ-1)，以减小近地面网格的偏斜角度，改善网格质量。在逐层向上调整优化的过程中，还可设置调整的阈值并使修正量逐渐减小，以避免出现负体积或网格尺寸变化过快；

以对网格编号为(i,j,k)的网格点的X坐标进行修正为例，Fortran语言程序示例为：

alph=atand( (Xc(i+1,j-1,k)-Xc(i,j-1,k))/(Yc(i+1,j-1,k)-Yc(i,j-1,k)) )

beta=atand( (Xc(i,j-1,k)-Xc(i-1,j-1,k))/(Yc(i,j-1,k)-Yc(i-1,j-1,k))*(-1.0) )

if(alph<0.0) then

alph=alph+180.0

endif

if(beta<0.0) then

beta=beta+180.0

endif

temp=(Yc(i,j,k)-Yc(i,j-1,k))*tand((alph-beta)/2.0)

Xc(i,j,k)=Xc(i,j-1,k)+temp

对X坐标的修正量与网格在垂向的高度以及特征角度(α-β)/2的正切值成正比，其中，α为程序中的alph，β为程序中的beta。特征角度指的是相邻的下层网格的相邻边的角平分线与Y轴的夹角，如图5所示，∠DBE即为特征角度，图中水平线段DE的长度即为对X坐标的原始修正量，实际修正量可以受到阈值等约束条件的限制。

正交修正前后的网格对比如图6和7所示（图6、7中左侧一列的数据为最大夹角，右侧一列的数据为对应夹角区间内的网格个数），经正交修正后，贴近地面的几层网格的最大偏斜角有明显的改善。

S36包括步骤：按照网格分辨率需求（即目标分辨率）对水平方向坐标数组Xc、Zc进行下采样，可采用间隔取点、插值等方式进行，获得CFD求解过程所需要的目标网格点的水平坐标X、Z；

例如，以X和Z方向分别每间隔istep和kstep个网格取点为例，Fortran语言程序示例为：

kk=-1

do k=0,NZ-1,kstep

kk=kk+1

do j=0,NY-1

ii=0

do i=0,NX-1,istep

X(1,ii,j,kk)=Xc(i,k)

Y(2,ii,j,kk)=Yc(i,j,k)

Z(3,ii,j,kk)=Zc(i,k)

ii=ii+1

enddo

S37包括步骤：将S35步骤中所采集到的目标网格点所对应的目标网格数组输出（优选为PLOT3D格式），其中，X和Z方向实际输出的网格点数为NXr和NZr，Y方向输出的网格点数为NY。

Fortran语言程序示例为：

open(10, file='plot3dgrid.xyz', form='unformatted')

write(10) NXr,NY,NZr

write(10) (((X(1,i,j,k),i=0,NXr-1),j=0,NY-1),k=0,NZr-1),&

(((Y(2,i,j,k),i=0,NXr-1),j=0,NY-1),k=0,NZr-1),&

(((Z(3,i,j,k),i=0,NXr-1),j=0,NY-1),k=0,NZr-1)

close(10)

网格划分程序生成的复杂地形表面网格及垂向网格分布如图8所示，复杂地形网格计算域整体视图如图9所示。

进一步地，在输出PLOT3D格式的目标网格数组之后，还可将PLOT3D格式转存为CGNS等通用的CFD求解格式。具体地，使用Pointwise等网格划分软件读取PLOT3D格式坐标数据，设置边界条件并导出为CGNS等通用的CFD求解格式。本发明方法借助在CFD领域应用广泛的Pointwise网格划分软件，可以直接读取HSMP程序（即上述网格划分程序）输出的PLOT3D格式文件，实现了便捷地将PLOT3D数据格式转存为CGNS等CFD通用格式，也可以根据CFD求解器类型自主选择需要导出的网格格式。

实施例二

如图10所示，本发明第二方面还提供了一种快速生成复杂地形表面高保真结构网格的***，包括：

数据获取模块10，被配置为用于获取具有地形数据的STL格式文件，所述地形数据包括：用于表示所述复杂地形表面的几何特征的多个三角形面片，所述三角形面片通过法线向量以及顶点坐标进行描述，其中，所述顶点坐标从投影坐标系中获取得到，所述投影坐标系包括：X方向和Z方向，以及与所述X方向、Z方向相互垂直的垂向Y方向；

在一些实施例中，本发明实施例中所选用的地形数据可以从GIS***导出， GIS是基于计算机技术对地理空间进行数据处理和数据应用的有用工具和技术。

当然，在另一些实施例中，也可以通过其他途径获取到地形表面的数据（如矢量地图数据），如通过open street map网站等。又或者，可以采用专业的GIS软件对高清地图图片进行矢量化从而得到对应地形的矢量地图数据。并通过数据格式转换方法（如转换软件等）将矢量地图数据导出为STL格式数据

数据提取模块20，被配置为用于采用坐标数值提取方法从所述地形数据中的目标范围内提取出多个坐标数据，并将所述多个坐标数据以离散点的形式进行存储，其中，所述坐标数据表示为（X，Y，Z）；

网格生成模块30，被配置为用于基于所述坐标数据生成具有规则排序的地面第一层网格数组，并从所述Y方向对所述地面第一层网格数组进行拉伸以创建形成目标网格数组；

其中，网格生成模块30包括：

二次排序单元，被配置为用于按照三维数组Zc中的坐标数值沿着Z方向的大小对所述三维数组Zc进行第二次排序，并同步变更三维数组Xc和Yc中对应坐标数值的位置，使得所述三维数组 Xc、Yc、Zc形成在Y方向上的地面第一层网格数组；

进一步地，在一些实施例中，该***还包括：导出模块，被配置为用于设定边界条件，将所述网格坐标点数组导出为用于模拟所述复杂地形表面的目标结构网格。

进一步地，在一些实施例中，网格生成模块30还包括：正交修正单元，被配置为用于在所述层网格中，从第二层网格开始，通过下层网格与相邻上层网格之间的夹角关系，分别在至少一个平面内对至少一个所述上层网格中的至少一个离散点的坐标进行修正。具体地修正方式可参见上述实施例。

本发明实施例中的***可以实现上述实施例一中的任一方法步骤，此处不再赘述。

可以理解的是，近地面区域的风力发电设备和超高层建筑（如特大型桥梁）的安全评估仅为本发明的一个优选应用场景。当然，本发明中的高保真结构网格生成方法也可以被广泛应用于其他风资源评估或其他相关领域，例如气象预报或者飞行安全分析等领域。

通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到上述实施例方法可借助软件加必需的通用硬件平台的方式来实现，当然也可以通过硬件，但很多情况下前者是更佳的实施方式。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质(如ROM/RAM、磁碟、光盘)中，包括若干指令用以使得一台计算机终端(可以是手机，计算机，服务器，或者网络设备等)执行本发明各个实施例所述的方法。

上面结合附图对本发明的实施例进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可做出很多形式，这些均属于本发明的保护之内。

Claims

1.一种快速生成复杂地形表面高保真结构网格的方法，其特征在于，包括步骤：

其中，S103包括步骤：

2.根据权利要求1所述的方法，其特征在于，在S35之前，还包括步骤：

；

所述离散点D修正后的坐标数值X'为：

；

所述离散点D修正后的坐标数值Z'为：

；

其中，和分别为ZY平面内的两个特征角∠CBD、∠ABD，离散点A₂的坐标为：，网格编号为，离散点B₂的坐标为：，网格编号为，离散点C₂的坐标为：，网格编号为，离散点D的坐标为：，网格编号为；

其中，所述网格编号用于表示所述坐标数据在二次排序之后的位置分布。

3.根据权利要求2所述的方法，其特征在于，在至少一个平面内对至少一个所述上层网格中的至少一个离散点的坐标进行修正之前，还包括步骤：

当所述总和小于预设的第一阈值时，则继续在当前的所述上层网格中进行修正，否则，停止在所述上层网格中进行修正；

和/或,

计算所述层网格中每一个已修正离散点的Z坐标值在Y方向上的修正量的总和；

4.根据权利要求2所述的方法，其特征在于，在至少一个平面内对至少一个所述上层网格中的至少一个离散点的坐标进行修正的步骤之前，还包括步骤：

当所述总和小于预设的第一层修正阈值时，则继续在当前的所述上层网格中进行修正，否则，停止在所述上层网格中进行修正；

和/或，

计算当前上层网格中每一个已修正的所述离散点的Z坐标值在Y方向上的修正值总和；

5.根据权利要求1所述的方法，其特征在于，还包括步骤：

6.根据权利要求1所述的方法，其特征在于，S31包括：

7.根据权利要求6所述的方法，其特征在于，NX为一维数组Xn中具有同一数值大小的坐标数值个数，NZ为沿Z方向的离散点的个数；其中，NZ=N/NX， N为所述有效坐标点的个数。

8.根据权利要求1所述的方法，其特征在于，S35包括：

9.根据权利要求1所述的方法，其特征在于，S36中采样的方法包括：邻点插值法，和/或间隔取点法。

10.一种快速生成复杂地形表面高保真结构网格的***，其特征在于，包括：

其中，网格生成模块包括：