CN115983163B - 基于sph算法的流域建模方法 - Google Patents

Abstract

本发明公开了一种基于SPH算法的流域建模方法，包括步骤S1：将流域划分成若干网格，并沿流域宽度方向将网格分为三个区域，包括靠近流域河床的两侧外部网格区域和位于流域中心的内部网格区域；步骤S2：流域由若干水分子微粒组成，建立水分子微粒受力模型，分析不同区域水分子微粒的受力情况，采用SPH算法得出水分子微粒的加速度关系式；步骤S3：每个网格中含有多个水分子微粒，给不同网格设置不同编号，建立网格哈希表，根据网格编号获取内部水分子微粒索引；然后根据水分子微粒的速度、加速度值，得出下一帧画面中水分子微粒所在位置，实现水流动态可视化模拟，效果更接近真实情况，为流域后续深入研究奠定良好的基础。

Description

基于SPH算法的流域建模方法

技术领域

本发明属于水文领域，具体涉及一种基于SPH算法的流域建模方法。

背景技术

水文监测是水文工作的重要组成部分，对我国水利规划、水工程建设管理、防汛抗早、水资源管理与保护具有重要意义。近年来，随着通信技术、人工智能技术等新技术的发展，以及全国中小河流水文监测系统工程、全国水文基础设施建设规划等重点项目的实施，我国水文水资源监测能力得到显著提升。

在对流域水文研究中，对流域进行模拟、分析也尤为重要。现有技术中，对流域的模拟研究中，大多数基于沿河流宽度方向水的流速相同这一基础，但是这样的设定与真实情况不符，导致模拟失真，进而影响后续的研究精准度。因为在真实环境中，由于流域两侧靠近河岸的的部分，水流受到河岸的阻力，导致流速较缓，而河流中间位置的流速更快。

发明内容

本发明提供了一种基于SPH算法的流域建模方法，解决现有技术中在对流域的模拟研究过程中，大多数研究基于沿河流宽度方向水的流速相同这一基础，但是这样的设定与真实情况不符，导致模拟失真的技术问题。

为了解决上述技术问题，本发明采用如下方案实现：

基于SPH算法的流域建模方法，包括如下步骤：

步骤S1：将流域划分成若干网格，并沿流域宽度方向将网格分为三个区域，包括靠近流域河床的两侧外部网格区域和位于流域中心的内部网格区域；

步骤S2：流域由若干水分子微粒组成，建立水分子微粒受力模型，分析不同区域水分子微粒的受力情况，采用SPH算法得出水分子微粒的加速度关系式；

步骤S3：每个网格中含有多个水分子微粒，给不同网格设置不同编号，建立网格哈希表，根据网格编号获取内部水分子微粒索引；然后根据水分子微粒的速度、加速度值，得出下一帧画面中水分子微粒所在位置，模拟流域中水的流动。

本申请中，基于水流受到河岸的阻力导致不同区域流速不同的真实状况，将流域沿其宽度方向将网格分为三个区域，分为靠近流域河床的两侧外部网格区域和位于流域中心的内部网格区域；对外部网格区域和内部网格区域的水流分别进行受力分析。根据不同区域水流的受力不同，得出不同区域水流下一时刻的速度、加速度，以及下一帧所在位置，实现水流动态可视化模拟，效果更接近真实情况，为流域后续深入研究奠定良好的基础。

进一步优化，所述步骤S2中，建立单个水分子微粒受力模型，并进行受力分析：

1)、靠近流域河床的两侧外部网格区域中的水分子微粒受力如下:

2)、流域中心的内部网格区域中的水分子微粒受力如下:

其中，

为水分子微粒的重力，g为重力加速度，ρ为密度；

为水分子微粒的受到的压力，/>

为水分子微粒从上一个位置流动到下一位置过程中所受压力之差；

为水分子微粒受到的粘性力，由水分子微粒之间的速度差产生，μ为粘度系数，/>

水分子微粒从初始位置流动到下一位置的流速差；

为水分子微粒受到的河床的阻力，由水流与河岸接触面摩擦产生；

其中，

λ表示水流与河岸的摩擦系数，n为与水流接触的河岸的粗糙度，e_t边界层厚度；S表示接触面积，Δz表示水分子微粒到河道接触面之间的垂直距离。

进一步优化，所述步骤S2中，采用SPH算法得出水分子微粒的加速度关系式，具体包括如下步骤：

步骤S2.1：定义光滑核函数，在流域中随机选择一个水分子微粒D，其半径为r；以水分子微粒D为中心，在光滑核半径h范围内有N个其他水分子微粒，半径依次为r₁,r₂,……，r_N，受周围其他N个水分子微粒作用后，则水分子微粒D的某项属性A为：

其中，j∈[1，N]，且为正整数，A_j为要累加的某种属性，m_j、ρ_j和

分别是周围其他水分子微粒的质量、密度和位置向量，/>

为水分子微粒D的位置向量，h为光滑核半径；W为光滑核函数；

步骤S2.2：分析单个水分子微粒与光滑核半径h范围内的其他水分子微粒之间的关系，具体为：

步骤S2.2.1：用密度ρ代替公式(3)中的A，设定所有水分子微粒的质量相同都为m，在三维空间中，

处的密度计算公式最终为：

依次用压力P、速度

代替公式(3)中的A，并与式(4)联合，得出光滑核半径h范围内水分子微粒所受重力、压力和粘性力；

步骤S2.2.2：计算靠近流域河床的两侧外部网格区域中水分子微粒所受阻力：

靠近流域河床的两侧外部网格区域中，沿流域宽度方向包含多层网格，不同层网格中的水分子微粒受到河床的阻力不同；具体为从河岸向河流中心，水分子微粒受到河床的阻力逐渐衰减，设定衰减系数β；

式中，k为水分子微粒所在边界层网格的索引，其中边界层最靠近河床壁面的那层网格索引为0，e_t为边界层厚度；

得出阻力产生的加速度：

步骤S2.2.2：计算出水分子微粒的加速度：

其中，流域中心的内部网格区域中水分子微粒的加速度：

靠近流域河床的两侧外部网格区域中的水分子微粒加速度：

进一步优化，所述步骤S3中，具体包括如下步骤：

步骤S3.1：首先创建单个水分子微粒的信息结构体，保存其位置、密度、压力、流体阻力、加速度、速度以及最终速度信息；其次，建立缓冲缓存类，保存所有的水分子微粒信息；并建立一个链接，将网格中的多个水分子微粒通过链接方式连接；

步骤S3.2：遍历链接，根据每个水分子微粒的受力，从而得出水分子微粒下一时刻的速度、加速度，以及下一帧画面中水分子微粒所在位置；伴随时间的推移，更新缓存数据，实现基于SPH的水分子微粒动态可视化；所有水分子微粒动态移动形成流域中水的流动。

进一步优化，靠近流域河床的两侧外部网格区域中，沿流域宽度方向包含多层网格，设定与河床接触的最外层网格编号为k＝0，次外层编号为1，依次递加；然后使用哈希表的方式，给每个空间网格赋予编号，根据网格编号在哈希表中获取水分子微粒索引，具体为：

首先建立三个网格哈希表，分别记为第一哈希表、第二哈希表和第三哈希表，其中，第一哈希表储存全部的网格信息，第二哈希表储存靠近流域河床的两侧外部网格区域中的网格信息，第三哈希表储存流域中心的内部网格区域的网格信息；

然后再根据靠近流域河床的两侧外部网格区域中的分层情况，设置一定数量的子哈希表，将外部网格区域中各层网格信息存储到对应编号的子哈希表中。

进一步优化，建立冲突链表，哈希冲突时，将网格中的多个水分子微粒通过链表方式连接。

进一步优化，所述网格边长设置为光滑核半径的2倍。

与现有技术相比，本发明具有如下有益效果：

本申请中，基于水流受到河岸的阻力导致不同区域流速不同的真实状况，将流域沿其宽度方向分为三个区域，分为靠近流域河床的两侧外部网格区域和位于流域中心的内部网格区域；对外部网格区域和内部网格区域的水流分别进行受力分析，根据不同区域水流的受力不同，得出不同区域水流下一时刻的速度、加速度，以及下一帧所在位置，实现水流动态可视化模拟，效果更接近真实情况，为流域后续深入研究奠定良好的基础。

附图说明

为了更清楚地说明本发明具体实施方式或现有技术中的技术方案，下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为发明所述基于SPH算法的流域建模方法的流程图；

图2为发明实施例中所模拟的戴村坝流域局部的初始位置模拟图；

图3为发明实施例中所模拟的戴村坝流域运行5秒后的局部模拟图。

具体实施方式

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

实施例一：

本实施例以戴村坝流域作为目标流域，戴村坝流域位于山东省东平县境内，为大汶河支流。如图1所示，基于SPH算法的流域建模方法，包括如下步骤：

步骤S1：将流域划分成若干网格，并沿流域宽度方向将网格分为三个区域，包括靠近流域河床的两侧外部网格区域和位于流域中心的内部网格区域；

步骤S2：流域由若干水分子微粒组成，建立水分子微粒受力模型，分析不同区域水分子微粒的受力情况，采用SPH算法得出水分子微粒的加速度关系式。SPH算法的基本设想，就是将连续的流体想象成一个个相互作用的水分子微粒，这些水分子微粒相互影响，共同形成了复杂的流体运动。每个水分子微粒拥有自己的特定属性，如位置，半径，所受重力黏力等等。因此在下一帧画面运动的时候，检索每个粒子周围一定范围内的水分子微粒，根据它们的位置(即计算出周围的微粒和当前微粒的间隔)重新计算黏力、压力和水流阻力，以及重力，根据受力计算新的加速度和速度，然后每一帧重复上述计算。因此在出给所有水分子微粒赋值初始位置及属性后，水分子微粒就能自动的开始动画直到最后平稳下来。

步骤S3：每个网格中含有多个水分子微粒，给不同网格设置不同编号，建立网格哈希表，根据网格编号获取内部水分子微粒索引；然后根据水分子微粒的速度、加速度值，得出下一帧画面中水分子微粒所在位置，模拟流域中水的流动。

本申请中，基于水流受到河岸的阻力导致不同区域流速不同的真实状况，将流域沿其宽度方向将网格分为三个区域，分为靠近流域河床的两侧外部网格区域和位于流域中心的内部网格区域；对外部网格区域和内部网格区域的水流分别进行受力分析。根据不同区域水流的受力不同，得出不同区域水流下一时刻的速度、加速度，以及下一帧所在位置，实现水流动态可视化模拟，效果更接近真实情况，为流域后续深入研究奠定良好的基础。

在本实施例中，所述步骤S2中，建立单个水分子微粒受力模型，并进行受力分析：

1)、靠近流域河床的两侧外部网格区域中的水分子微粒受力如下:

2)、流域中心的内部网格区域中的水分子微粒受力如下:

其中，

为水分子微粒的重力，g为重力加速度，ρ为密度；

为水分子微粒的受到的压力，/>

为水分子微粒从上一个位置流动到下一位置过程中所受压力之差。由于存在压力差水才会源源不断的流动，压力差大小等于压力场的梯度，方向由压力高的区域指向压力低的区域。

为水分子微粒受到的粘性力，由水分子微粒之间的速度差产生，μ为粘度系数，/>

水分子微粒从初始位置流动到下一位置的流速差。流体内部，快速流动的部分会施加类似于剪切力的作用力到速度慢的部分，这个力的大小跟流体的粘度系数μ以及速度差有关。

为水分子微粒受到的河床的阻力，由水流与河岸接触面摩擦产生。

其中，

λ表示水与河岸的摩擦系数，n为与水流接触的河岸的粗糙度，e_t边界层厚度；S表示接触面积，Δz表示水分子微粒到河道接触面之间的垂直距离。

在本实施例中，所述步骤S2中，采用SPH算法得出水分子微粒的加速度关系式，具体包括如下步骤：

步骤S2.1：定义光滑核函数，在流域中随机选择一个水分子微粒D，其半径为r；以水分子微粒D为中心，在光滑核半径h范围内有N个其他水分子微粒，半径依次为r₁,r₂,……，r_N，受周围其他N个水分子微粒作用后，则水分子微粒D的某项属性A为：

其中，j∈[1，N]，且为正整数，A_j为要累加的某种属性，m_j、ρ_j和

分别是周围其他水分子微粒的质量、密度和位置向量，/>

为水分子微粒D的位置向量，h为光滑核半径；W为光滑核函数。

步骤S2.2：分析单个水分子微粒与光滑核半径h范围内周围的其他水分子微粒之间的关系，具体为：

步骤S2.2.1：用密度ρ代替公式(3)中的A，设定所有水分子微粒的质量相同都为m，在三维空间中，

处的密度计算公式最终为：

依次用压力P、速度

代替公式(3)中的A，并与式(4)联合，得出光滑核半径h范围内水分子微粒所受重力、压力和粘性力。

步骤S2.2.2：计算靠近流域河床的两侧外部网格区域中水分子微粒所受阻力：

靠近流域河床的两侧外部网格区域中，沿流域宽度方向包含多层网格，不同层网格中的水分子微粒受到河床的阻力不同。具体为从河岸向河道中心，水分子微粒受到河床的阻力逐渐衰减，设定衰减系数β；

式中，k为水分子微粒所在边界层网格的索引，其中边界层最靠近河床壁面的那层网格索引为0，e_t为边界层厚度。

在SPH算法里，流体的质量由流体单元的密度决定，所以一般用密度代替水分子微粒的质量，即F^resistance＝ρ_ia^resistance；

则有：

由于河岸静止，所有速度差记为水粒子的速度：Δu＝u_i；

边界层中每一个子层离河岸的平均距离为：Δz＝(2k-1)h；

由上式联立得出阻力产生的加速度：

步骤S2.2.2：计算出水分子微粒的加速度：

其中，流域中心的内部网格区域中的水分子微粒加速度：

骤S2.2.2：计算出水分子微粒的加速度：

其中，流域中心的内部网格区域中水分子微粒的加速度：

靠近流域河床的两侧外部网格区域中的水分子微粒加速度：

在本实施例中，所述步骤S3中，具体包括如下步骤：

步骤S3.1：首先创建单个水分子微粒的信息结构体，保存其位置、密度、压力、流体阻力、加速度、速度以及最终速度信息；其次，建立缓冲缓存类，保存所有的水分子微粒信息；并建立一个链接，将网格中的多个水分子微粒通过链接方式连接。

步骤S3.2：遍历链接，根据每个水分子微粒的受力，从而得出水分子微粒下一时刻的速度、加速度，以及下一帧画面中水分子微粒所在位置；伴随时间的推移，更新缓存数据，实现基于SPH的水分子微粒动态可视化；所有水分子微粒动态移动形成流域中水的流动。

在本实施例中，水分子微粒受到重力、压力、粘性力，河道阻力，且会与河岸产生碰撞，运动一段时间后粒子位置会发生变化，图2为初始位置模拟图；图3为5秒后的模拟图。

图2中靠近河岸的A点位置、速度变化如表1所示，图2中靠近河流中心的B点位置、速度变化如表2所示。

表1靠近河岸的A点位置、速度变化数据

表2靠近河流中心的B点位置、速度变化数据

在本实施例中，由于计算每个水分子微粒的状态时，都需要获得在它光滑核半径内(邻域内)的所有水分子微粒信息。如果遍历每个水分子微粒，计算欧式距离的话，开销过于庞大，因此本申请中流域将空间划分成多个均匀的网格，从而在获取每个水分子微粒的邻域粒子时，只需要访问该水分子微粒的周围网格获取邻域粒子即可。靠近流域河床的两侧外部网格区域中，沿流域宽度方向包含多层网格，设定与河床接触的最外层网格编号为k＝0，次外层编号为1，依次递加；然后使用哈希表的方式，给每个空间网格赋予编号，根据网格编号在哈希表中获取水分子微粒索引，具体为：先建立三个网格哈希表，分别记为第一哈希表、第二哈希表和第三哈希表，其中，第一哈希表储存全部的网格信息，第二哈希表储存靠近流域河床的两侧外部网格区域中的网格信息，第三哈希表储存流域中心的内部网格区域的网格信息；然后再根据靠近流域河床的两侧外部网格区域中的分层情况，设置一定数量的子哈希表，将外部网格区域中各层网格信息存储到对应编号的子哈希表中。建立冲突链表，哈希冲突时，将网格中的多个水分子微粒通过链表方式连接。

流体系统需要一个水分子微粒缓存类，来保存所有的水分子微粒信息。空间网格划分类，包含空间网格的编号以及哈希表信息。由于在每帧迭时，所有的水分子微粒信息都会改变，每个水分子微粒信息的邻域信息都会改变，所以还需要建立一个邻接表类，用来保存每个水分子微粒信息的邻域粒子的信息。最后将上述SPH算法中流域相关内容整合到一个流体系统类当中。

SPH算法包括初始化函数_init，密度和压力计算函数_computerPressure，加速度计算函数_computerForce和移动粒子函数_advance。模拟的水流动画在tick帧函数中进行。计算每帧水分子微粒位置后，更新posData缓存数据，则实现了基于SPH的粒子动态可视化。

在本实施例中，所述网格边长设置为光滑核半径的2倍。这样在二维空间中最多只遍历4个网格就可以获取邻域水分子微粒。在三维空间中，则最多只需要遍历2*2*2，8个网格，就可以获取邻域水分子微粒。具体为：当网格点在三维空间的顶点，则遍历1个网格；网格点在三维空间的线上，则遍历2个网格；网格点在三维空间的面上，则遍历4个网格；网格点在三维空间内部，则需要遍历8个网格。

以上述依据本发明的理想实施例为启示，通过上述的说明内容，相关工作人员完全可以在不偏离本项发明技术思想的范围内，进行多样的变更以及修改。本项发明的保护范围并不局限于说明书上的内容，必须要根据权利要求范围来确定其保护范围。

Claims (6)

1.基于SPH算法的流域建模方法，其特征在于，包括如下步骤：

步骤S1：将流域划分成若干网格，并沿流域宽度方向将网格分为三个区域，包括靠近流域河床的两侧外部网格区域和位于流域中心的内部网格区域；

步骤S2：流域由若干水分子微粒组成，建立水分子微粒受力模型，分析不同区域水分子微粒的受力情况，采用SPH算法得出水分子微粒的加速度关系式；

建立单个水分子微粒受力模型，并进行受力分析，具体为：

1)、靠近流域河床的两侧外部网格区域中的水分子微粒受力如下:

2)、流域中心的内部网格区域中的水分子微粒受力如下:

其中，

为水分子微粒的重力，g为重力加速度，ρ为密度；

为水分子微粒的受到的压力，/>

为水分子微粒从上一个位置流动到下一位置过程中所受压力之差；

为水分子微粒受到的粘性力，由水分子微粒之间的速度差产生，μ为粘度系数，/>

水分子微粒从初始位置流动到下一位置的流速差；

为水分子微粒受到的河床的阻力，由水流与河岸接触面摩擦产生；

其中，

λ表示水流与河岸的摩擦系数，n为与水流接触的河岸的粗糙度，e_t边界层厚度；S表示接触面积，Δz表示水分子微粒到河道接触面之间的垂直距离；

步骤S3：每个网格中含有多个水分子微粒，给不同网格设置不同编号，建立网格哈希表，根据网格编号获取内部水分子微粒索引；然后根据水分子微粒的速度、加速度值，得出下一帧画面中水分子微粒所在位置，模拟流域中水的流动。

2.根据权利要求1所述的基于SPH算法的流域建模方法，其特征在于，所述步骤S2中，采用SPH算法得出水分子微粒的加速度关系式，具体包括如下步骤：

步骤S2.1：定义光滑核函数，在流域中随机选择一个水分子微粒D，其半径为r；以水分子微粒D为中心，在光滑核半径h范围内有N个其他水分子微粒，半径依次为r₁,r₂,……，r_N，受周围其他N个水分子微粒作用后，则水分子微粒D的某项属性A为：

其中，j∈[1，N]，且为正整数，A_j为要累加的某种属性，m_j、ρ_j和

分别是周围其他水分子微粒的质量、密度和位置向量，/>

为水分子微粒D的位置向量，h为光滑核半径；W为光滑核函数；

步骤S2.2：分析单个水分子微粒与光滑核半径h范围内的其他水分子微粒之间的关系，具体为：

步骤S2.2.1：用密度ρ代替公式(3)中的A，设定所有水分子微粒的质量相同都为m，在三维空间中，

处的密度计算公式最终为：

依次用压力P、速度

代替公式(3)中的A，并与式(4)联合，得出光滑核半径h范围内水分子微粒所受重力、压力和粘性力；

步骤S2.2.2：计算靠近流域河床的两侧外部网格区域中水分子微粒所受阻力：

靠近流域河床的两侧外部网格区域中，沿流域宽度方向包含多层网格，不同层网格中的水分子微粒受到河床的阻力不同；具体为从河岸向河流中心，水分子微粒受到河床的阻力逐渐衰减，设定衰减系数β；

式中，k为水分子微粒所在边界层网格的索引，其中边界层最靠近河床壁面的那层网格索引为0，e_t为边界层厚度；

得出阻力产生的加速度：

步骤S2.2.2：计算出水分子微粒的加速度：

其中，流域中心的内部网格区域中水分子微粒的加速度：

靠近流域河床的两侧外部网格区域中的水分子微粒加速度：

3.根据权利要求2所述的基于SPH算法的流域建模方法，其特征在于，所述步骤S3中，具体包括如下步骤：

步骤S3.1：首先创建单个水分子微粒的信息结构体，保存其位置、密度、压力、流体阻力、加速度、速度以及最终速度信息；其次，建立缓冲缓存类，保存所有的水分子微粒信息；并建立一个链接，将网格中的多个水分子微粒通过链接方式连接；

步骤S3.2：遍历链接，根据每个水分子微粒的受力，从而得出水分子微粒下一时刻的速度、加速度，以及下一帧画面中水分子微粒所在位置；伴随时间的推移，更新缓存数据，实现基于SPH的水分子微粒动态可视化；所有水分子微粒动态移动形成流域中水的流动。

4.根据权利要求3所述的基于SPH算法的流域建模方法，其特征在于，靠近流域河床的两侧外部网格区域中，沿流域宽度方向包含多层网格，设定与河床接触的最外层网格编号为k＝0，次外层编号为1，依次递加；然后使用哈希表的方式，给每个空间网格赋予编号，根据网格编号在哈希表中获取水分子微粒索引，具体为：

首先建立三个网格哈希表，分别记为第一哈希表、第二哈希表和第三哈希表，其中，第一哈希表储存全部的网格信息，第二哈希表储存靠近流域河床的两侧外部网格区域中的网格信息，第三哈希表储存流域中心的内部网格区域的网格信息；

然后再根据靠近流域河床的两侧外部网格区域中的分层情况，设置一定数量的子哈希表，将外部网格区域中各层网格信息存储到对应编号的子哈希表中。

5.根据权利要求4所述的基于SPH算法的流域建模方法，其特征在于，建立冲突链表，哈希冲突时，将网格中的多个水分子微粒通过链表方式连接。

6.根据权利要求2所述的基于SPH算法的流域建模方法，其特征在于，所述网格边长设置为光滑核半径的2倍。

Images