CN109102558A - 一种基于非结构三角网的水流流线快速生成方法 - Google Patents

Abstract

本发明涉及一种基于非结构三角网的水流流线快速生成方法，包括如下步骤：1）非结构三角网拓扑结构预处理，构建“点‑单元”拓扑关系表；2）设置水流流线搜索起始点，并找出起始点所在单元编号；3）采用变搜索步长+水波型外扩的递归搜索方式，向起始点的上下游方向进行水流流线追踪；4）按上游水流流线搜索终点—起始点—下游水流流线搜索终点的顺序对水流流线点进行排序。本发明在进行水流流线生成时，采用了水波型外扩的递归搜索方式，最大程度地减小了在进行水流流线搜索时进行坐标点与网格单元的关系判断的次数，显著地提高了水流流线的搜索效率；同时采用了动态搜索步长的方法，在网格密集区域自动减小搜索步长确保了网格密集区域水流流线生成的精度，在网格稀疏区域自动加大搜索步长避免了大量无效计算，提高了水流流线的生成效率。

Description

一种基于非结构三角网的水流流线快速生成方法

技术领域

本发明涉及一种基于非结构三角网的水流流线快速生成方法，属于计算机应用技术领域。

背景技术

水流流线被定义为流场中每一点上都与速度矢量相切的曲线，是表述矢量场模型的一种有效方式。采用水流流线进行矢量场表述的方式既弥补了箭头表达缺失的连续性，又补充了颜色映射无法表征方向的缺陷，被广泛运用于各类水利项目的流场分析，是水利数值模拟后处理的重要功能之一。

在实际运用中，水流流线通常是在水动力学数值模型模拟得到的各离散单元的流速、流向信息的基础上，采用特定的水流流线追踪算法由用户指定的起始点向其上下游方向追踪得到。模拟区域离散方法不同，其对应的水流流线追踪算法也不尽相同。近年来，随着计算机性能的不断提高和有限元、有限体积算法的发展，几何适应性较强的非结构三角网逐渐取代了结构网格，越来越多地用于水动力数值模拟中模拟区域的离散。

非结构三角网的一个显著特点是单元间不存在明确的拓扑规则，即无法直接从拓扑关系中找出两个单元间的关系。另外，由于非结构三角网的几何特性，同一套网格内的网格尺寸通常不统一，特别是在实际运用过程中经常出现需要对重点关注区域进行局部网格加密的情况，同一套网格内的网格尺寸通常会出现一个甚至多个量级上的差别。针对非结构三角网的这些特点，目前的水流流线生成算法中存在以下两个主要的技术问题：

（1）由于非结构三角网单元间不存在明确的拓扑规则，在进行水流流线搜索时需要耗费大量的计算资源用于遍历所有单元来进行搜索点与网格单元位置关系的判断。特别是在网格单元数量大、水流流线生成量多时，进行水流流线生成的耗时长、效率低；

（2）针对非结构三角网在实际运用过程中，同一套网格中网格尺寸通常差异较大的特点，传统的固定步长水流流线搜索方法有着较大的局限性，搜索步长取值太小会增加大量的计算量影响搜索效率，搜索步长取值太大则会直接影响水流流线的搜索精度，尤其是在重点关注的网格加密区。固定搜索步长的取值难以在计算量和计算精度之间进行平衡。

鉴于以上基于非结构三角网的水流流线生成中存在的问题，寻找一种在保证精度的条件下，能显著提高基于非结构三角网的水流流线生成速度的方法，对于水利数值模拟的实际运用具有重要的意义。

发明内容

本发的目的是提出一种基于非结构三角网的水流流线生成方法，在不牺牲搜索精度的前提下，解决在大网格单元数量的非结构三角网基础上进行水流流线搜索速度慢的问题。

本发明为解决上述技术问题，提供一种基于非结构三角网的水流流线快速生成方法。具体包括以下步骤：

步骤一：非结构三角网拓扑结构预处理，构建“点-单元”拓扑关系表。

通常，非结构三角网拓扑关系文件包括：①点列表，包含点的编号及对应坐标信息；②边列表，包含边的编号、构成边的点编号以及边的类型；③单元列表，包含单元的编号、构成单元的点编号。在这些基础拓扑关系的基础上，为提高水流流线搜索效率，需要在进行开始搜索水流流线前，对基础非结构三角网的拓扑关系列表进行扩展，增加“点-单元”拓扑关系列表，即非结构三角网各点与所有包含该点的单元的关系列表。

步骤二：设置水流流线搜索起始点，并找出起始点所在单元编号。

水流流线搜索需要一个搜索起始点，起始点可以由用户输入坐标或在屏幕上点选。起始点选择后，遍历非结构三角网所有单元，采用重心法对起始点坐标与各单元位置关系进行判断，找到起始点所在的单元编号，并将起始点标记为当前搜索点。

采用重心法进行点与单元位置关系的判断方法如下：

；

；

式中：P为要进行判断的点，A、B、C分别为三角形单元的三个顶点，为A点到B点的向量，为A点到C点的向量，为A点到P点的向量；

当1≥a≥0、1≥b≥0、a+b≤1三个条件同时满足时，P在三角形ABC内部（含位于边界上），否则P在三角形ABC外部。

步骤三：采用变搜索步长+水波型外扩的递归搜索方式，依次向起始点的上游和下游方向进行水流流线追踪。

起始点设置完成后将该点作为构成水流流线的一个点，依次向该点的上游和下游方向分别进行水流流线的追踪，直至上下游方向均遇到当前点不在网格范围内或当前点所在单元流速为0的情况，搜索结束，按上游水流流线搜索终点—起始点—下游水流流线搜索终点的顺序对搜索所得的各点进行排序。

水流流线的追踪的具体步骤如下：

步骤31，为保证水流流线搜索精度，采用插值的方法提取当前点X、Y方向流速信息。根据当前点所在单元三个顶点的X、Y方向流速与当前点距离所在单元三个顶点的距离，插值出当前点的X、Y方向流速u、v。具体公式如下：

；；

式中，u₁、u₂、u₃分布为单元三个顶点的X方向流速；v₁、v₂、v₃分别为单元三个顶点的Y方向流速；d₁、d₂、d₃为当前点与所在单元三个顶点的距离；u、v为当前点的X、Y方向流速。

计算得到流速u、v值后，对u、v值进行判断：

①若u、v值其中至少有一个不等于0，即当前点流速不为0，则记录当前点的坐标信息，然后进入步骤32；

②若u、v值均等于0，即当前点流速为0，表示当前搜索方向的水流流线到达终点。此时，若当前搜索方向为上游方向，则重复步骤三，由起始点开始向下游方向开始追踪水流流线。若当前搜索方向为下游方向，表示上下游方向水流流线搜索均已到达终点，则结束步骤三，进入步骤四。

步骤32，根据当前点的X、Y方向流速u、v，计算出当前搜索方向（上游/下游）的单位向量。

下游方向：；

上游方向：；

式中：为当前搜索方向的单位向量；u、v为当前点的X、Y方向流速；

步骤33，根据单元尺寸设置初始搜索步长。根据当前搜索点所在单元节点坐标，计算当前单元的平均边长，取单元的平均边长的1/2作为搜索步长Δs的初始值：

；

式中：Δs为当前搜索步长；为所在单元的平均边长。

步骤34，根据当前搜索点的X、Y坐标、X、Y方向的流速u、v及搜索步长Δs计算出下一个点的坐标，将计算出的点位作为当前搜索点。下一个点的坐标计算公式为：

（X_new，Y_new）=（X，Y）+Δs*；

式中，X_new、Y_new为下一个点的坐标；X、Y为当前点的坐标；Δs为当前搜索步长；为当前搜索方向的单位向量。

步骤35，以上一个点所在单元为起始搜索范围，采用水波型外扩的递归算法搜索下一个点所在单元。

创建一个以搜索点坐标和搜索网格范围列表为参数的搜索函数F和一个记录当前递归调用层数的变量C。搜索函数F的作用是寻找搜索点所在的网格单元，具体执行方式为：

①当搜索点落在搜索网格范围内的单元时，返回所在单元编号；

②当递归调用层数大于5层时，则认为当前搜索点已经超出非结构三角网的范围，当前搜索方向的水流流线搜索到达终点，此时，若当前搜索方向为上游方向，则重复步骤三，由起始点开始向下游方向开始追踪水流流线。若当前搜索方向为下游方向，表示上下游方向水流流线搜索均已到达终点，则结束步骤三，进入步骤四。

③搜索点不在当前搜索网格范围内且当前递归调用层数小于等于5层时，则利用步骤一构建的“点-单元”拓扑关系表，将位于当前搜索网格范围外侧且与当前搜索网格范围内单元有共同顶点的全部单元作为新的搜索网格范围，以新的搜索网格范围列表作为输入参数进行搜索函数F的递归调用，进行下一层搜索。搜索网格范围水波型外扩方式如图2所示，图中1号网格为起始搜索网格，1、2、3标识的网格分别表示第1~3层搜索网格范围。

步骤36，若步骤35的搜索函数返回的是单元编号，则利用本次搜索的递归调用层数对搜索步长进行动态调整，然后进行下一个点搜索。

本专利中，搜索的递归调用层数反映的是以上一个水流流线点所在网格单元为中心，搜索到下一个点需要外扩的网格单元层数，为了平衡搜索精度和搜索效率，外扩网格层数取值为2较为合理，即搜索得到的两个相邻水流流线点所在的网格单元相邻。因此，搜索步长Δs根据以下规则进行调整：

；

式中：Δs为调整前搜索步长；为调整后搜索步长；C为本次搜索的递归调用层数。

调整搜索步长后，重复步骤31~步骤35直至起始点上下游方向的水流流线搜索均到达终点。

步骤四：对水流流线点进行排序。

由于本方法中水流流线搜索方向是由起始点向上下游方向分别搜索，因此在水流流线搜索完成后，需要对搜索得到的水流流线点按上游水流流线搜索终点—起始点—下游水流流线搜索终点的顺序进行排序，得到最终水流流线。

本发明相对于现有技术具有如下的优点：

（1）本发明中的水流流线生成方法在进行水流流线生成时，采用了水波型外扩的递归搜索方式，最大程度地减小了在进行水流流线搜索时进行坐标点与网格单元的关系判断的次数，显著地提高了水流流线的搜索效率；

（2）本发明中的水流流线生成方法在进行水流流线生成时，采用了动态搜索步长的方法，在网格密集区域自动减小搜索步长确保了网格密集区域水流流线生成的精度，在网格稀疏区域自动加大搜索步长避免了大量无效计算，提高了水流流线的生成效率。

附图说明

图1 为本发明的方法流程图。

图2为搜索网格范围水波型外扩示意图。

图3 为实施例中用于水流流线生成示例的非结构三角网。

图4为 图3矩形区域放大之后的非结构三角网。

图5 为实施例最终生成的水流流线。

具体实施方式

下面结合实施例对本发明内容做进一步的说明：图1 为本发明的方法流程图。

本发明的一个实施例，以一个基于非结构三角网的二维水动力模型模拟结果为基础进行水流流线生成，用于验证本发明的可行性和有效性。

本实施例选取的用于水流流线生成的河网模型实例的建模河网对象位于广州市番禺区，北纬22°52'37.1"〜22°57'15"，东经113°15'41"〜113°29'47"，由紫坭河水道、龙湾河、市桥水道、沙湾水道构成的一段河网。其中，图3为河网模拟范围的非结构三角网剖分，图4为图3中矩形区域的局部放大图。

步骤一：非结构三角网拓扑结构预处理，构建“点-单元”拓扑关系表。

对二维水动力模型的非结构三角网进行处理，遍历非结构三角网的点列表和单元列表，生成“点-单元”拓扑关系列表，如下所示：

表1 非结构三角网的“点-单元”拓扑关系列表

点编号	单元编号
		1	1、2、3、4、6、8
2	6、8、7、5、10、9
		……	……
25768	49867、49868、49861、49862、48796、48795

步骤二：设置水流流线搜索起始点，并找出起始点所在单元编号。

根据分析需要，由用户输入水流流线起始点X=38443030，Y=2534730；遍历非结构三角网所有单元，采用重心法判断得该点所在单元编号为31040，将该点设置为当前搜索点。

步骤三：采用变搜索步长+水波型外扩的递归搜索方式，依次向起始点的上游和下游方向进行水流流线追踪。

首先提取当前点所在单元的三个顶点的坐标及X、Y方向流速u₁、u₂、u₃和v₁、v₂、v₃。根据坐标信息计算当前点与三个顶点的距离d₁、d₂、d₃。以起始点为例，提取及计算结果如下表所示：

表2 当前点所在单元信息提取结果

点编号	X，Y坐标	与当前点距离	u	v
					24649	38443057，2534730	27.00	0.58	-0.03
24558	38443027，2534732	3.61	0.59	-0.02
					24645	38443041，2534711	21.95	0.57	-0.03

距离及流速信息提取后，采用以下公式对当前点的速度信息进行插值：

； ；

计算的当前点的u=0.586，v=0.22。

对u、v值进行判断：①若u、v值其中至少有一个不等于0，即当前点流速不为0，则记录当前点的坐标信息，然后进入步骤32；②若u、v值均等于0，即当前点流速为0，表示当前搜索方向的水流流线到达终点。此时，若当前搜索方向为上游方向，则重复步骤三，由起始点开始向下游方向开始追踪水流流线。若当前搜索方向为下游方向，表示上下游方向水流流线搜索均已到达终点，则结束步骤三，进入步骤四。根据计算出的u、v值，属于前一类情况，记录当前点坐标信息。

下一步，根据当前搜索点的X、Y方向流速u、v，根据以下公式计算出当前搜索方向（上游/下游）方向的单位向量：

下游方向：；

上游方向：；

此处以计算上游方向单位向量为例，计算得到上游方向单位向量为（-0.999，0.038）。

下一步，根据当前搜索点所在的单元节点坐标，计算出当前单元的平均边长，取单元的平均边长的1/2作为搜索步长Δs的初始值。经过计算，当前单元的平均边长为26.71，初始搜索步长Δs=26.71/2=13.36。

下一步，根据当前搜索点的X、Y坐标、X、Y方向的流速u、v及搜索步长Δs，利用以下公式计算出下一个点的坐标。

（X_new，Y_new）=（X，Y）+Δs*；

计算得到下一个点的坐标为（38443016.65，2534730.51），将计算出的点位作为当前搜索点。

接下来，采用水波型外扩的递归算法搜索下一个点所在单元，以上一个点所在31040单元为起始搜索范围。首先判断当前搜索点（38443016.65，2534730.51）是否落在31040单元内，判断结果为否且当前递归调用层数小于5，则检索步骤一生成的“点-单元”拓扑关系表，将位于31040单元外侧且与31040单元有共同顶点的全部单元作为新的搜索网格范围。遍历新网格搜索范围内单元，依次判断当前搜索点是否落在当前单元内，判断结果为当前搜索点落在30926单元内，返回所在单元编号。本次搜索过程两层递归搜索的网格编号见表3。

表3 搜索网格范围

递归调用层数	搜索网格范围包含的单元编号	当前点是否落在该范围内
			1	31040	否
2	30925、31039、31161、30926、31035、31036、31156、31155、31162、30921、31279、31285	是，落在30926单元内

最后根据当前递归调用层数，利用以下公式对搜索步长进行调整：

；

由于本次递归调用层数C=2,计算得到Δs_new=Δs=13.36。搜索步长调整后，重复步骤三，直至上下游方向的水流流线搜索均到达终点。

步骤四：对水流流线点进行排序。

对搜索得到的水流流线点按上游水流流线搜索终点—起始点—下游水流流线搜索终点的顺序进行排序，得到如图5所示的最终水流流线。

本发明未进一步展开详细说明的内容均属于技术领域的公知常识，以上实施案例是本发明技术方案的一个实例，本发明的公开范围并不局限于此，任何对于本发明的技术方案的等同替换，都落在本发明的保护范围之内。

Claims (2)

1.一种基于非结构三角网的水流流线快速生成方法，其特征在于，包含如下步骤：

步骤一：非结构三角网拓扑结构预处理，构建“点-单元”拓扑关系表；

步骤二：设置水流流线搜索起始点，并找出起始点所在单元编号；

步骤三：采用变搜索步长与水波型外扩的递归搜索方式，依次向起始点的上游和下游方向进行水流流线追踪；

步骤四：按上游水流流线搜索终点-搜索起始点-下游水流流线搜索终点的顺序对水流流线点进行排序。

2.根据权利要求1所述的基于非结构三角网的水流流线快速生成方法，其特征在于，所述步骤三中的递归搜索方式包括：

步骤31：根据当前点所在单元三个顶点的X、Y方向流速与当前点距离所在单元三个顶点的距离，采用插值的方法插值出当前点的X、Y方向流速u、v；计算公式如下：

； ；

式中，u1、u2、u3分布为单元三个顶点的X方向流速；v1、v2、v3分别为单元三个顶点的Y方向流速；d1、d2、d3为当前点与所在单元三个顶点的距离；u、v为当前点的X、Y方向流速；

计算得到流速u、v值后，对u、v值进行判断：

①u、v值其中至少有一个不等于0，则记录当前点的坐标信息，执行步骤32；

②若u、v值均等于0，则进一步根据当前搜索方向进行判断：如果为上游方向，则重复执行步骤三；如果为下游方向，则执行步骤四；

步骤32：根据当前点的X、Y方向流速u、v，计算出当前搜索方向上游、下游的单位向量：

下游方向：；

上游方向：；式中：为当前搜索方向的单位向量；u、v为当前点的X、Y方向流速；

步骤33：根据单元尺寸设置初始搜索步长：根据当前搜索点所在单元节点坐标，计算当前单元的平均边长，取单元的平均边长的1/2作为搜索步长Δs的初始值：；式中：Δs为当前搜索步长；为所在单元的平均边长；

步骤34：根据当前搜索点的X、Y坐标、X、Y方向的流速u、v及搜索步长Δs计算出下一个点的坐标，将计算出的点位作为当前搜索点；下一个点的坐标计算公式为：（X_new，Y_new）=（X，Y）+Δs*；式中，X_new、Y_new为下一个点的坐标；X、Y为当前点的坐标；Δs为当前搜索步长；为当前搜索方向的单位向量；

步骤35：以上一个点所在单元为起始搜索范围，采用水波型外扩的递归算法搜索下一个点所在单元：包括创建一个以搜索点坐标和搜索网格范围列表为参数的搜索函数F和一个记录当前递归调用层数的变量C，具体方式为：

①当搜索点落在搜索网格范围内的单元时，返回所在单元编号；

②当递归调用层数大于5层时，则认为当前搜索点已经超出非结构三角网的范围，当前搜索方向的水流流线搜索到达终点，此时，若当前搜索方向为上游方向，则重复步骤三，由起始点开始向下游方向开始追踪水流流线；若当前点搜索方向为下游方向，表示上下游方向水流流线搜索均已到达终点，则进入步骤四；

③搜索点不在当前搜索网格范围内且当前递归调用层数小于等于5层时，则利用步骤一构建的“点-单元”拓扑关系表，将位于当前搜索网格范围外侧且与当前搜索网格范围内单元有共同顶点的全部单元作为新的搜索网格范围，以新的搜索网格范围列表作为输入参数进行搜索函数F的递归调用，进行下一层搜索；

步骤36：若步骤35的搜索函数返回的是单元编号，则利用本次搜索的递归调用层数对搜索步长进行动态调整，返回执行步骤31；

搜索步长Δs根据以下规则进行调整：；式中Δs为调整前搜索步长；为调整后搜索步长；C为本次搜索的递归调用层数。