CN109711525A - 一种用于sph算法的邻近粒子搜索方法及系统 - Google Patents

Abstract

本发明公开了一种用于SPH算法的邻近粒子搜索方法及系统，涉及计算流体力学技术领域。本发明通过构建预置单胞网格并计算每个粒子所在位置的网格单元编号；将粒子映射到网格单元中，获取每个粒子所对应的网格单元编号，根据网格单元编号找出粒子的相邻网格单元；逐个比较每个粒子与相邻网格单元中的粒子的距离并建立粒子对。该方法可以有效提高软件的鲁棒性。通过算例验证，使用本发明进行流体力学问题的数值计算仿真是可行的，有效提高了计算效率。

Description

一种用于SPH算法的邻近粒子搜索方法及系统

技术领域

本发明涉及计算流体力学技术领域，具体是涉及一种用于SPH算法的邻近粒子搜索方法及系统。

背景技术

随着计算机技术和CFD计算方法的迅速发展，SPH法逐渐用于流体数值计算过程中，但是传统的SPH法计算规模大，计算效率低，因此，在现有有效研究手段缺乏的情况下，开展SPH法改进算法研究工作来提高计算效率显得尤为重要。

对于流体域内任一个粒子，其核函数有一个支持域，在支持域内的所有粒子都会与该粒子发生相互作用，都是该粒子的临近粒子。临近粒子会随计算时间的变化而改变，因此在每一个迭代步都必须进行临近粒子的搜索算法，才能确定并计算粒子对之间的相互作用。相对于整个计算域来说，单个粒子的支持域非常小，因此其临近粒子数相对于整个计算域粒子数来说也很少，因SPH方法是无网格法，粒子之间无连接关系，因此如何高效地从数量庞大、散乱无规则的粒子群中找出它们的临近粒子，是亟待解决的重要问题。

粒子的排列顺序对计算速度也有较明显的影响，因为在计算粒子对之间相互作用的时候，需要访问两个粒子的数据。CPU访问数据时，总线往往会将数据所在的内存及附近的数据调入缓存区，若两粒子在内存中的位置相距较远，则CPU从缓存区取得数据的几率下降。对于整个计算域来说，若每个粒子与它的配对粒子相隔较远，则会耗费大量的CPU处理时间，导致程序的整体计算效率降低。

发明内容

本发明的目的是为了克服上述背景技术的不足，提供一种用于SPH算法的邻近粒子搜索方法及系统，在大规模无网格流体计算时，提高计算效率。

本发明提供一种用于SPH算法的邻近粒子搜索方法，包括以下步骤：

建立所有粒子区域的笛卡尔坐标系，构建预置单胞网格并计算每个粒子所在位置的网格单元编号；

将粒子映射到网格单元中，获取每个粒子所对应的网格单元编号，根据网格单元编号找出粒子的相邻网格单元；

逐个比较每个粒子与相邻网格单元中的粒子的距离并建立粒子对。

在上述方案的基础上，所述构建预置单胞网格并计算每个粒子所在位置的网格单元编号，具体包括以下步骤：

遍历所有粒子，计算出粒子所在坐标系的最大坐标值(X_max,Y_max,Z_max)、最小坐标值(X_min,Y_min,Z_min)和粒子的最大光滑长度h_max；

计算预置单胞网格的边长d，ε为容差，计算公式为：

d＝2h_max+ε (1)

将粒子所在坐标系用边长d等分成大小相同的网格单元，计算网格单元编号最大值：(i_max,j_max,k_max)：

计算网格单元总数Cell Num：

Cell Num＝i_max*j_max*k_max (3)；

计算每个粒子所在位置的网格单元编号(i,j,k)，将粒子对

号入座：

在上述方案的基础上，所述将粒子映射到网格单元中，获取每个粒子所对应的网格单元编号，根据网格单元编号找出粒子的相邻网格单元，具体包括以下步骤：

将所有粒子根据计算出的每个粒子的网格单元编号映射到网格单元中；

获取每个粒子所对应的网格单元编号(i,j,k)；

根据网格单元位置，找出粒子的相邻网格单元，2维情况下相邻网格单元有8个，3维情况下相邻网格单元有26个。

在上述方案的基础上，所述逐个比较每个粒子与相邻网格单元中的粒子的距离并建立粒子对，具体包括以下步骤：

逐个比较与相邻网格单元中的粒子的距离，若距离小于2倍光滑长度，则与其建立粒子对，若距离大于2倍光滑长度，则认为两个粒子之间没有相互作用，不建立粒子对。

在上述方案的基础上，该方法还包括以下步骤，逐个比较每个粒子与相邻网格单元中的粒子的距离并建立粒子对之前，对粒子在内存中的存储顺序进行排序，排序方法为：

对每个粒子，计算：

根据(i_p,j_p,k_p)，按位置顺序计算每个粒子的Morton key，转换成二进制码后，按照编码插空顺序排列；

将所有粒子按照Morton key值从小到大调整在内存中的存储顺序，最后按照Z-curve order算法重新排序。

在上述方案的基础上，当第一次粒子排序之后，后续的SPH计算中再次进行粒子排序时，比较当前粒子排序结果和上一次粒子排序结果，仅对顺序发生改变的粒子进行内存数据移动操作。

本发明还提供一种用于SPH算法的邻近粒子搜索系统，包括：

网格编号计算模块，其用于：建立所有粒子区域的笛卡尔坐标系，构建预置单胞网格并计算每个粒子所在位置的网格单元编号；

相邻网格获取模块，其用于：将粒子映射到网格单元中，获取每个粒子所对应的网格单元编号，根据网格单元编号找出粒子的相邻网格单元；

粒子对建立模块，其用于：逐个比较每个粒子与相邻网格单元中的粒子的距离并建立粒子对。

在上述方案的基础上，该系统还包括粒子排序模块，其用于：逐个比较每个粒子与相邻网格单元中的粒子的距离并建立粒子对之前，对粒子在内存中的存储顺序进行排序，排序方法为：

根据(i_p,j_p,k_p)，按位置顺序计算每个粒子的Morton key，转换成二进制码后，按照编码插空顺序排列；

将所有粒子按照Morton key值从小到大调整在内存中的存储顺序，最后按照Z-curve order算法重新排序。

与现有技术相比，本发明的优点如下：

(1)本发明通过构建预置单胞网格并计算每个粒子所在位置的网格单元编号；将粒子映射到网格单元中，获取每个粒子所对应的网格单元编号，根据网格单元编号找出粒子的相邻网格单元；逐个比较每个粒子与相邻网格单元中的粒子的距离并建立粒子对。该方法可以有效提高软件的鲁棒性。通过算例验证，使用本发明进行流体力学问题的数值计算仿真是可行的，有效提高了计算效率。

(2)本发明逐个比较每个粒子与相邻网格单元中的粒子的距离并建立粒子对之前，根据Morton key对粒子在内存中的存储顺序进行排序。采用本发明的方法进行粒子排序后，粒子与临近粒子的内存位置相距较近，可明显改观程序的计算效率，对并行计算速度的影响尤其显著。

(3)当第一次粒子排序之后，后续的SPH计算中再次进行粒子排序时，比较当前粒子排序结果和上一次粒子排序结果，仅对顺序发生改变的粒子进行内存数据移动操作。而不必对所有的粒子都重新操作一次，大大提升计算效率。

附图说明

图1是本发明实施例的预置单胞网格和粒子配对算法示意图；

图2是本发明实施例的计算Morton Key的方法示意图；

图3是本发明实施例的粒子排序前的示意图；

图4是本发明实施例的粒子排序后的示意图；

图5是本发明实施例的不同水粒子尺寸；

图6是本发明实施例的水粒子在0.01s时刻的半球运动状态；

图7是本发明实施例的水粒子在0.02s时刻的半球运动状态；

图8是本发明实施例的水粒子在0.03s时刻的半球运动状态。

具体实施方式

下面结合附图及具体实施例对本发明作进一步的详细描述。

实施例1：

本发明实施例提供一种用于SPH算法的邻近粒子搜索方法，包括以下步骤：

参见图1所示，建立所有粒子区域的笛卡尔坐标系，构建预置单胞网格并计算每个粒子所在位置的网格单元编号；

将粒子映射到网格单元中，获取每个粒子所对应的网格单元编号，根据网格单元编号找出粒子的相邻网格单元；

逐个比较每个粒子与相邻网格单元中的粒子的距离并建立粒子对。

所述构建预置单胞网格并计算每个粒子所在位置的网格单元编号，具体包括以下步骤：

遍历所有粒子，计算出粒子所在坐标系的最大坐标值(X_max,Y_max,Z_max)、最小坐标值(X_min,Y_min,Z_min)和粒子的最大光滑长度h_max；

计算预置单胞网格的边长d，ε为容差，计算公式为：

d＝2h_max+ε (1)

将粒子所在坐标系用边长d等分成大小相同的网格单元，计算网格单元编号最大值：(i_max,j_max,k_max)：

计算网格单元总数Cell Num：

Cell Num＝i_max*j_max*k_max (3)；

计算每个粒子所在位置的网格单元编号(i,j,k)，将粒子对

号入座：

所述将粒子映射到网格单元中，获取每个粒子所对应的网格单元编号，根据网格单元编号找出粒子的相邻网格单元，具体包括以下步骤：

将所有粒子根据计算出的每个粒子的网格单元编号映射到网格单元中；

获取每个粒子所对应的网格单元编号(i,j,k)；

根据网格单元位置，找出粒子的相邻网格单元，2维情况下相邻网格单元有8个，3维情况下相邻网格单元有26个。

所述逐个比较每个粒子与相邻网格单元中的粒子的距离并建立粒子对，具体包括以下步骤：

逐个比较与相邻网格单元中的粒子的距离，若距离小于2倍光滑长度，则与其建立粒子对，若距离大于2倍光滑长度，则认为两个粒子之间没有相互作用，不建立粒子对。

该方法可以有效提高软件的鲁棒性。通过算例验证，使用本发明进行流体力学问题的数值计算仿真是可行的，有效提高了计算效率。

作为优选的实施方式，该方法还包括以下步骤，逐个比较每个粒子与相邻网格单元中的粒子的距离并建立粒子对之前，对粒子在内存中的存储顺序进行排序，参见图2所示，本发明实施例中和了Morton key(摩顿秘钥)和Z-curve order(Z线性排序)的方法对粒子排序。排序方法为：

根据(i_p,j_p,k_p)，按位置顺序计算每个粒子的Morton key，转换成二进制码后，按照编码插空顺序排列；

将所有粒子按照Morton key值从小到大调整在内存中的存储顺序，最后按照Z-curve order算法重新排序。

本发明实施例逐个比较每个粒子与相邻网格单元中的粒子的距离并建立粒子对之前，根据Morton key对粒子在内存中的存储顺序进行排序。图3为排序前水粒子的顺序，图4为排序后水粒子的顺序，排序前其粒子之间的内存间距在局部一些位置较大，采用该方法进行排序后，可以明显看出，从左下角开始，粒子顺序呈Z形分布，粒子与临近粒子的内存位置相距较近。采用本发明的方法进行粒子排序后，粒子与临近粒子的内存位置相距较近，可明显改观程序的计算效率，对并行计算速度的影响尤其显著。

作为优选的实施方式，当第一次粒子排序之后，后续的SPH计算中再次进行粒子排序时，比较当前粒子排序结果和上一次粒子排序结果，仅对顺序发生改变的粒子进行内存数据移动操作。而不必对所有的粒子都重新操作一次，大大提升计算效率。

实施例2：

本发明实施例提供一种用于SPH算法的邻近粒子搜索系统，包括：

网格编号计算模块，其用于：建立所有粒子区域的笛卡尔坐标系，构建预置单胞网格并计算每个粒子所在位置的网格单元编号；

相邻网格获取模块，其用于：将粒子映射到网格单元中，获取每个粒子所对应的网格单元编号，根据网格单元编号找出粒子的相邻网格单元；

粒子对建立模块，其用于逐个比较每个粒子与相邻网格单元中的粒子的距离并建立粒子对。

该系统可以有效提高软件的鲁棒性。通过算例验证，使用本发明进行流体力学问题的数值计算仿真是可行的，有效提高了计算效率。

作为优选的实施方式，该系统还包括粒子排序模块，其用于：逐个比较每个粒子与相邻网格单元中的粒子的距离并建立粒子对之前，对粒子在内存中的存储顺序进行排序，排序方法为：

根据(i_p,j_p,k_p)，按位置顺序计算每个粒子的Morton key，转换成二进制码后，按照编码插空顺序排列；

将所有粒子按照Morton key值从小到大调整在内存中的存储顺序，最后按照Z-curve order算法重新排序。

本发明实施例逐个比较每个粒子与相邻网格单元中的粒子的距离并建立粒子对之前，根据Morton key对粒子在内存中的存储顺序进行排序。采用本发明的方法进行粒子排序后，粒子与临近粒子的内存位置相距较近，可明显改观程序的计算效率，对并行计算速度的影响尤其显著。

实施例3：

参见图5～图8所示，验证模型采用圆柱形的水域，冲击物为半球形，进行自由落体着水试验。对多核并行计算方法进行了测试，还对临近粒子搜索和排序算法进行了测试。

当调用单个核进行求解，并且不使用粒子排序算法的情况下，0.01m水粒子网格的仿真计算用时80分钟41秒。当调用4个核进行求解，并且不使用粒子排序算法的情况下，0.01m水粒子网格的仿真计算用时37分钟23秒。可见，使用并行计算方法可以大幅度提高SPH法的计算效率。

当调用单个核进行求解，并且使用粒子排序算法的情况下，0.01m水粒子网格的仿真计算用时58分钟18秒。因此，使用粒子排序算法也可以有效提高SPH法的计算效率。通过该算例验证，证明了使用本发明的算法进行该类问题仿真是可行的，结果是可信的。

本领域的技术人员可以对本发明实施例进行各种修改和变型，倘若这些修改和变型在本发明权利要求及其等同技术的范围之内，则这些修改和变型也在本发明的保护范围之内。说明书中未详细描述的内容为本领域技术人员公知的现有技术。

Claims (8)

1.一种用于SPH算法的邻近粒子搜索方法，其特征在于，包括以下步骤：

建立所有粒子区域的笛卡尔坐标系，构建预置单胞网格并计算每个粒子所在位置的网格单元编号；

将粒子映射到网格单元中，获取每个粒子所对应的网格单元编号，根据网格单元编号找出粒子的相邻网格单元；

逐个比较每个粒子与相邻网格单元中的粒子的距离并建立粒子对。

2.如权利要求1所述的方法，其特征在于：所述构建预置单胞网格并计算每个粒子所在位置的网格单元编号，具体包括以下步骤：

遍历所有粒子，计算出粒子所在坐标系的最大坐标值(X_max,Y_max,Z_max)、最小坐标值(X_min,Y_min,Z_min)和粒子的最大光滑长度h_max；

计算预置单胞网格的边长d，ε为容差，计算公式为：

d＝2h_max+ε (1)

将粒子所在坐标系用边长d等分成大小相同的网格单元，计算网格单元编号最大值：(i_max,j_max,k_max)：

计算网格单元总数Cell Num：

Cell Num＝i_max*j_max*k_max (3)

计算每个粒子所在位置的网格单元编号(i,j,k)，将粒子对号入座：

i＝int[(X-X_min)/d]+1 (4)。

j＝int[(Y-Y_min)/d]+1

k＝int[(Z-Z_min)/d]+1

3.如权利要求1所述的方法，其特征在于：所述将粒子映射到网格单元中，获取每个粒子所对应的网格单元编号，根据网格单元编号找出粒子的相邻网格单元，具体包括以下步骤：

将所有粒子根据计算出的每个粒子的网格单元编号映射到网格单元中；

获取每个粒子所对应的网格单元编号(i,j,k)；

根据网格单元位置，找出粒子的相邻网格单元，2维情况下相邻网格单元有8个，3维情况下相邻网格单元有26个。

4.如权利要求1所述的方法，其特征在于：所述逐个比较每个粒子与相邻网格单元中的粒子的距离并建立粒子对，具体包括以下步骤：

逐个比较与相邻网格单元中的粒子的距离，若距离小于2倍光滑长度，则与其建立粒子对，若距离大于2倍光滑长度，则认为两个粒子之间没有相互作用，不建立粒子对。

5.如权利要求1所述的方法，其特征在于：该方法还包括以下步骤，逐个比较每个粒子与相邻网格单元中的粒子的距离并建立粒子对之前，对粒子在内存中的存储顺序进行排序，排序方法为：

对每个粒子，计算：

根据(i_p,j_p,k_p)，按位置顺序计算每个粒子的Morton key，转换成二进制码后，按照编码插空顺序排列；

将所有粒子按照Morton key值从小到大调整在内存中的存储顺序，最后按照Z-curveorder算法重新排序。

6.如权利要求1所述的方法，其特征在于：当第一次粒子排序之后，后续的SPH计算中再次进行粒子排序时，比较当前粒子排序结果和上一次粒子排序结果，仅对顺序发生改变的粒子进行内存数据移动操作。

7.一种用于SPH算法的邻近粒子搜索系统，其特征在于，包括：

网格编号计算模块，其用于：建立所有粒子区域的笛卡尔坐标系，构建预置单胞网格并计算每个粒子所在位置的网格单元编号；

相邻网格获取模块，其用于：将粒子映射到网格单元中，获取每个粒子所对应的网格单元编号，根据网格单元编号找出粒子的相邻网格单元；

粒子对建立模块，其用于：逐个比较每个粒子与相邻网格单元中的粒子的距离并建立粒子对。

8.如权利要求7所述的系统，其特征在于：该系统还包括粒子排序模块，其用于：逐个比较每个粒子与相邻网格单元中的粒子的距离并建立粒子对之前，对粒子在内存中的存储顺序进行排序，排序方法为：

根据(i_p,j_p,k_p)，按位置顺序计算每个粒子的Morton key，转换成二进制码后，按照编码插空顺序排列；

将所有粒子按照Morton key值从小到大调整在内存中的存储顺序，最后按照Z-curveorder算法重新排序。