CN102930087B - 一种模拟仿真技术中的相邻粒子搜索方法 - Google Patents

Abstract

本发明公开了一种模拟仿真技术中的相邻粒子快速搜索方法。本发明方法对传统的PIB搜索法进行了重大改进，克服了PIB搜索法对粒子点分布的依赖，具有比PIB搜索法和树形搜索法更高的计算效率；该方法适用于各种商用计算机仿真系统，如：ANSYS、LS-DYNA、ABAQUS等，可嵌入到计算机仿真系统中以提高相邻粒子搜索的计算效率，从而提高软件模拟大规模工程问题的能力。

Description

一种模拟仿真技术中的相邻粒子搜索方法

技术领域

本发明涉及工程技术中的数值仿真领域，尤其涉及一种模拟仿真技术中的相邻粒子快速搜索方法。

背景技术

汽车碰撞、飞机降落以及爆炸冲击等工程实际问题的模拟仿真都存在计算量巨大的问题，其中接触界面的接触搜索、粒子法（如：光滑粒子流体动力学方法、离散元法、物质点法等）的相邻粒子搜索耗时较多是导致模拟仿真计算量大的主要原因之一。资料调研显示，目前已有的相邻粒子搜索方法包括：直接搜索法、链表搜索法、树形搜索法、PIB（Point-In-Box）搜索法等，现有技术中的参考文献1中对树形搜索法进行了公开，而参考文献2对PIB搜索法进行了公开。其中，直接搜索法效率太低，只对简单问题适用；链表搜索法在定区域的情况下效率很高，但在变区域的情况下效率低；树形搜索法是一种非常稳定、高效的算法；PIB搜索法的效率较高，而且非常节省内存，但其搜索效率对粒子分布很敏感。但以高效的树形搜索法应用于光滑粒子流体动力学方法的相邻粒子搜索为例，其搜索耗时仍然占整个仿真计算时间的50%以上。因此，必须发明一种更高效的搜索方法，以提高计算机模拟仿真的计算效率，使模拟仿真更好的解决工程实际问题。

所有已公开的关于计算机仿真方面的专利，包括：200710092856.6一种塑性成形数值模拟方法、200710171864.X管道中含有不可凝结气体的蒸汽冷凝的数值模拟方法、200810008558.9数值模拟结果显示程序、方法和系统、200610026093.0板柱结构中板柱节点的数值模拟方法、200710051543.6复杂渗控结构的渗流问题的SVA数值模拟方法、200710040000.4冲压模具结构分析数值模拟方法、200710035112.0基于数值模拟的高强铝、镁合金等温挤压方法、201110120406.X船舶与海洋工程防腐系统数值模拟及优化方法、201110088396.6一种级配碎石三轴试验的数值模拟方法、201110131761.7数值模拟金属薄板成形中预测颈缩破坏的方法和系统、201110194080.5一种铸件宏观偏析数值模拟的方法、201110226317.3选择性进流水温平抑装置及其水温数值模拟预报方法、201110231013.6粉沙质和淤泥质海岸泥沙运动数值模拟方法、201110327595.8木材复杂各向异性本构关系模型的数值模拟方法、201010290825.3一种基于数值模拟的汽车覆盖件回弹误差补偿方法、201110223231.5高强度钢板温热成形数值模拟方法、200810023420.6主梁断面气动自激力的全过程数值模拟方法、200810036138.1点焊连接失效数值模拟系统等。虽然专利较多，但模拟仿真技术中的相邻粒子快速搜索方法方面的相关专利目前还未有公开。

本发明主要针对以往搜索方法进行相邻粒子搜索效率低的问题，创新性的提出了一种条形化PIB搜索方法，该方法可有效提高工程实际问题中的模拟仿真计算效率。

[参考文献]

[1]Liu G R,Liu M B.Smoothed particle hydrodynamics:a meshfree particle method.Singapore:World Scientific,2003.

[2]Swegle J W.Search Algorithm.Sandia National Laboratories,1992.

发明内容

本发明的目的在于解决模拟仿真中相邻粒子搜索计算效率低的问题，提出一种可由计算机仿真系统实施的相邻粒子快速搜索方法，从而提高计算机仿真系统模拟复杂问题时的计算效率。

根据本发明的一个方面，提供一种模拟仿真技术中的相邻粒子搜索方法，所述方法包括以下步骤：

步骤1：将点集划分为一系列条形的子集；

步骤2：分别对每个子集内的点排序；

步骤3：搜索给定盒子内的点；以及

步骤4：判断盒子内的粒子是否与对象粒子形成相邻粒子对，从而搜索到所有相邻粒子。

优选地，在所述步骤1中，将点集占据的最小区域在垂直于条形方向的方向上进行分割，以得到若干小的条形区域，每一个小的条形区域中的点即构成为一个子集，其中条形方向是点集宽度最大的方向，即坐标跨度最大的方向。

优选地，在所述步骤2中，对每个子集中的点按各个坐标方向的坐标值排序，点排序的结果是索引数组和序号数组，其中该索引数组按局部排序序号的升序存储各个点集中点的局部编号，该序号数组按照局部编号的升序存储各个点集中点的局部排序序号，其中局部编号指点在子集中的编号，局部排序序号指点在子集中的排序序号。

优选地，在所述步骤3中，先确定可能包含位于盒子内的点的子集，再利用PIB搜索法逐个确定这些子集中位于盒子内的点，最后将在每个子集中搜索到的点合并即得到盒子内的所有点。

优选地，在所述步骤1中，当三维空间的条形区域的尺寸为Δ_s、条形的方向为z方向时，将点集占据的区域在x和y方向上进行分割，其分割数分别为：

tns^x＝int[(x_max-x_min)/Δ_s]+1和tns^y＝int[(y_max-y_min)/Δ_s]+1       (1)

公式（1）中：x_max、x_min分别为点的最大和最小x坐标，y_max、y_min分别为点的最大和最小y坐标；

将点集占据的区域在x和y方向上分割后，得到tns＝tns^x×tns^y个小的条形区域，对每个条形区域按如下方式编号：

ns(ns^x,ns^y)=(ns^y-1)tns^x+ns^x          (2)

公式（2）中：ns为条形区域的编号，ns^x和ns^y分别为条形区域在x和y方向的序号；

对于任意给定的点i，该点所在的条形区域在x和y方向上的序号按下式计算：

ns_i ^x＝min(int[(x_i-x_min)/Δ^x _s]+1,tns^x)

                                     (3)    ns_i ^y＝min(int[(y_i-y_min)/Δ^y _s]+1,tns^y)

公式（3）中： ${\mathrm{\Δ}}_s^x=(x_{\max }-x_{\min })/{\mathrm{tns}}^x$ 和 ${\mathrm{\Δ}}_s^y=(y_{\max }-y_{\min })/{\mathrm{tns}}^y$ 分别为x和y方向的精确的条形尺寸；

将按公式(3)计算的序号代入公式(2)，即得到点i的所在的条形区域的编号ns_i，采用上述方法确定每一点所在的条形区域的编号，所有在同一个条形区域中的点构成一个子集，子集的编号为所在条形区域的编号。

优选地，在所述步骤1中，在划分子集时构造以下数组来记录子集的相关信息从而方便后续的搜索操作：长度为总点数N的数组Strip和Ndsort，其中Strip存储每一个点的所在的子集编号，Ndsort按点所在子集的编号的升序存储点的编号；长度为tns的数组Nstrip和Npoint，其中Nstrip记录每一个子集中的点数，Npoint记录每一个子集中的第一个点在Ndsort中的位置，其中，构造这些数组的流程如下：

(1)初始化数组：Nstrip=0；

(2)根据公式（2）确定每一个点i所在的子集编号ns_i；

(3)存储点i的子集编号到数组Strip：Strip(i)=ns_i；

(4)修改子集ns_i中的点数：Nstrip(ns_i)=Nstrip(ns_i)+1；

(5)计算每一个子集j中的第一个点在Ndsort中的位置：

Npoint(1)=1,Npoint(j)=Npoint(j-1)+Nstrip(j-1)；

(6)初始化数组：Nstrip=0；

(7)将点按其所在的子集的编号的升序存储到Ndsort中：Ndsort(Nstrip(Strip(i))+Npoint(Strip(i)))=i，Nstrip(Strip(i))=Nstrip(Strip(i))+1。

优选地，在所述步骤2中，对每个子集中点的排序的具体流程如下：

(1)初始化j=1；

(2)确定子集j中的第一个和最后一个点在Ndsort中的位置：

jsta=Npoint(j),jend=Npoint(j)+Nstrip(j)-1；

(3)形成子集j的坐标数组（lx，ly，lz）：

lx(1～Nstrip(j))=x(Ndsort(jsta～jend))，ly(1～Nstrip(j))=y(Ndsort(jsta～jend))，lz(1～Nstrip(j))=z(Ndsort(jsta～jend))；

其中，上面三个式子中的波浪线“～”表示在取值范围内逐一取值；

(4)对lx、ly和lz分别排序，形成局部索引数组（I′^x、I′^y、I′^z）和局部序号数组（R′^x、R′^y、R′^z），其中I′^x(i)存储局部排序序号为i的点的局部编号，R′^x(i)存储局部编号为i的点的局部排序序号，I′^y、I′^z和R′^y、R′^z的含义分别与I′^x和R′^x类似；

(5)分别将局部索引数组和局部序号数组转化为索引数组和序号数组：

I^x(jsta～jend)=I′^x(1～Nstrip(j))，I^y(jsta～jend)=I′^y(1～Nstrip(j))，

I^z(jsta～jend)=I′^z(1～Nstrip(j))，R^x(jsta～jend)=R′^x(1～Nstrip(j))，

R^y(jsta～jend)=R′^y(1～Nstrip(j))，R^z(jsta～jend)=R′^z(1～Nstrip(j))。

(6)判断j是否等于tns，如果不是则j=j+1，返回步骤（2）继续循环排序，如果是则结束排序。

优选地，在所述步骤3中，以对象粒子a为中心构造边长为4h_a的正方体盒子k，其中h_a是粒子a的光滑长度，

搜索给定盒子k内的点的具体流程如下：

(1)计算包含盒子k的x和y方向边界的条形区域在x和y方向上的序号：

上式中：x_lk和x_tk为盒子k的两个x方向边界的位置,y_lk和y_tk为两个y方向边界的位置；

(2)x和y方向的序号分别在Istrip_min～Istrip_max和Jstrip_min～Jstrip_max之间的条形区域对应的子集可能包含位于盒子k内的点，对这些子集中的点进行搜索；

(3)合并上述步骤(2)中在每一个子集中搜索到的点，即得到点集中位于给定盒子k内的所有点。

优选地，在所述步骤3的上述步骤（2）中，对这些子集按照以下方式搜索：

(2.1)初始化Istrip=Istrip_min，开始循环1；

(2.2)初始化Jstrip=Jstrip_min，开始循环2；

(2.3)根据公式（2）计算x和y方向的序号分别为Istrip和Jstrip的条形区域对应的子集的编号ns；

(2.4)确定子集ns中第一个和最后一个点的位置：

ista=Npoint(ns)，iend=Npoint(ns)+Nstrip(ns)-1；

(2.5)形成局部索引数组和排序数组：

I′^x(1～Nstrip(ns))=I^x(ista～iend)，I′^y(1～Nstrip(ns))=I^y(ista～iend)，

I′^z(1～Nstrip(ns))=I^z(ista～iend)，R′^x(1～Nstrip(ns))=R^x(ista～iend)，

R′^y(1～Nstrip(ns))=R^y(ista～iend)，R′^z(1～Nstrip(ns))=R^z(ista～iend)；

(2.6)构造点列表：对局部索引数组采用二分法进行搜索，得到子集ns中在各个方向上的第一个和最后一个位于盒子内的点的局部排序序号（和）；局部排序序号在之间的点构成x方向的点列表，同样根据和分别确定y和z方向的点列表；

(2.7)寻找点列表的交集：选择点数最少的点列表，对其中的点逐个检查，判断其在其他方向的局部排序序号是否位于相应方向的第一个和最后一个点的排序序号之间：当z方向的点列表中的点数最少，则逐个检查局部编号为的点；如果且则局部编号为i′的点在盒子k内，否则不在其中；得到点的局部编号i′后，其全局编号由数组Ndsort得到，即i＝Ndsort(i′+Npoint(ns)-1)；

(2.8)判断Jstrip是否等于Jstrip_max，如果不是则Jstrip=Jstrip+1，返回步骤（2.3）继续循环2，如果是则结束循环2，转入步骤(2.9)；

(2.9)判断Istrip是否等于Istrip_max，如果不是则Istrip=Istrip+1，返回步骤（2.2）继续循环1，如果是则结束循环1从而结束步骤（2）的搜索。

优选地，在所述步骤4中，以如下方式来判断盒子内的粒子是否与对象粒子a形成相邻粒子对：

如果搜索到的粒子b满足：

a<b，d_ab<2h_a， $d_{\mathrm{ab}}<2{\stackrel{\&OverBar;}{h}}_{\mathrm{ab}}---\left(6\right)$

或者

a>b，d_ab≥2h_b， $d_{\mathrm{ab}}<2{\stackrel{\&OverBar;}{h}}_{\mathrm{ab}}---\left(7\right)$

则粒子a和粒子b记录为相邻粒子对，其中a和b均为粒子编号，上式中d_ab指粒子a、b间的距离，其中h_a是粒子a的光滑长度，h_b是粒子b的光滑长度。

从上述技术方案可见，本发明提出的相邻粒子快速搜索方法对传统的PIB搜索法进行了重大改进，克服了PIB搜索法对粒子点分布的依赖，具有比PIB搜索法和树形搜索法更高的计算效率；该方法适用于各种商用计算机仿真系统，如：ANSYS、LS-DYNA、ABAQUS等，可嵌入到计算机仿真系统中以提高相邻粒子搜索的计算效率，从而提高软件模拟大规模工程问题的能力。

附图说明

下面将通过参照附图详细描述本发明的优选实施例，使本领域的普通技术人员更清楚本发明的上述及其它特征和优点，附图中：

图1是本发明中的条形化PIB搜索法在二维空间中的实施示意图；

图2是本发明中索引数组和序号数组的存储示意图；

图3是二维空间中对象粒子的相邻粒子示意图；

图4是二维空间中条形分布的点集示意图；

图5是二维空间中方形分布的点集示意图；

图6是二维空间中不同点数和点的分布情况下的搜索耗时示意图；

图7是三维空间中条形分布的点集示意图；

图8是三维空间中方形分布的点集示意图；

图9是三维空间中不同点数和点的分布情况下的搜索耗时示意图；

图10是三维空间中不同点数情况下PIB搜索法与条形化PIB搜索法的搜索耗时比示意图；

图11是本发明中的条形化PIB搜索法与树形搜索法搜索相邻粒子的耗时示意图；

图12是本发明中的条形化PIB搜索法与树形搜索法搜索相邻粒子的耗时比示意图；

图13是卵形钢弹斜侵彻铝板的计算模型示意图；

图14是应用不同搜索法时一个光滑粒子法的计算循环的耗时比较示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，以下举实施例对本发明进一步详细说明。

为了提高PIB搜索法的性能，减小对粒子点分布形式的敏感性，利用PIB搜索法对条形分布点集搜索效率高和稳定的特点，本发明提出了一种条形化PIB搜索法，并将其应用于相邻粒子搜索。需要说明的是，在本发明中，“点”和“粒子”的含义是相同的，仅仅是本领域技术人员在不同处理过程中对粒子的习惯叫法而已。

本发明方法的基本思想是：先把点集划分为一系列条形的子集，分别对每个子集内的点排序，当搜索给定盒子（本发明中盒子即在粒子的周围构建一个边界盒，边界盒的边为直线，并且与体系的坐标系平行）内的点时，先确定可能包含位于盒子内的点的子集，再利用PIB搜索法逐个确定这些子集中位于盒子内的点，最后将在每个子集中搜索到的点合并即得到盒子内的所有点（粒子），之后就可以判断盒子内的粒子是否与对象粒子形成相邻粒子对，从而搜索到所有相邻粒子。

本发明该方法的具体步骤如下：

步骤1：将点集划分条形子集

本步骤用于将点集划分为一系列条形的子集。

划分条形子集首先需要确定条形区域的尺寸和方向。条形区域的尺寸对最终的搜索效率有影响，后文中在对实施例1的介绍中说明了条形尺寸Δ_s的确定方法。对于条形区域的方向，本发明中方法选择点集宽度最大的方向（即坐标跨度最大的方向），这样有利于计算效率的提高。划分子集的方法是：将点集占据的最小区域（即包含点集的最小长方体（对于三维空间而言）或矩形（对于二维空间而言））在垂直于条形方向的方向上进行分割，得到若干小的条形区域，每一个小的条形区域中的点即构成为一个子集。为示意方便，图1示出了本发明中的条形化PIB搜索法在二维空间中的实施示意图，而下面的具体实例则对条形化PIB搜索法在三维空间中的实施情况进行了说明。

假设条形的尺寸为Δ_s、条形的方向为z方向，则将点集占据的区域在x和y方向上进行分割，其分割数分别为：

tns^x＝int[(x_max-x_min)/Δs]+1和tns^y＝int[(y_max-y_min)/Δ_s]+1         (1)

以上公式（1）中：x_max、x_min分别为点的最大和最小x坐标，y_max、y_min分别为点的最大和最小y坐标。

将点集占据的区域在x和y方向上分割后，可以得到tns＝tns^x×tns^y个小的条形区域，对每个条形区域按如下方式编号：

ns(ns^x,ns^y)=(ns^y-1)tns^x+ns^x              (2)

以上公式（2）中：ns为条形区域的编号，ns^x和ns^y分别为条形区域在x和y方向的序号。

对于任意给定的点i，该点所在的条形区域在x和y方向上的序号按以下公式（3）计算：

ns_i ^x＝min(int[(x_i-x_min)/Δ^x _s]+1,tns^x)

                              (3)    ns_i ^y＝min(int[(y_i-y_min)/Δ^y _s]+1,tns^y)

以上公式（3）中： ${\mathrm{\&Delta;}}_s^x=(x_{\max }-x_{\min })/{\mathrm{tns}}^x$ 和 ${\mathrm{\&Delta;}}_s^y=(y_{\max }-y_{\min })/{\mathrm{tns}}^y$ 分别为x和y方向的精确的条形尺寸。

将按公式(3)计算的序号代入公式(2)，即得到点i的所在的条形区域的编号ns_i。采用上述方法，确定每一点所在的条形区域的编号，所有在同一个条形区域中的点构成一个子集，子集的编号为所在条形区域的编号。

为方便后续的搜索，在划分子集时需要构造以下数组来记录子集的相关信息：长度为N（总的点数）的数组Strip和Ndsort，其中Strip存储每一个点的所在的子集编号，Ndsort按点所在子集的编号的升序存储点的编号；长度为tns的数组Nstrip和Npoint，其中Nstrip记录每一个子集中的点数，Npoint记录每一个子集中的第一个点在Ndsort中的位置。构造这些数组的流程如下：

(1)初始化数组：Nstrip＝O；

(2)根据公式(2)确定每一个点i所在的子集编号ns_i；

(3)存储点i的子集编号到数组Strip：Strp(i)＝ns_i；

(4)修改子集ns_i中的点数：Nstrip(ns_i)＝Nstrip(ns_i)+1；

(5)计算每一个子集j中的第一个点在Ndsort中的位置：

Npoint(1)＝1，Npoint(j)＝Npoint(j-1)+Nstrip(j-1)；

(6)初始化数组：Nstrip＝O；

(7)将点按其所在的子集的编号的升序存储到Ndsort中：

Ndsort(Nstrip(Strip(i))+Npoint(Strip(i)))＝i，Nstrip(Str(i))＝Nstrip(Strip(i))+1。

步骤2：分别对每个子集内的点排序

划分完子集后，对每个子集中的点按各个坐标方向的坐标值排序。点排序的结果是索引数组(I^x、I^y、I^z)和序号数组(R^x、R^y、R^z)。图2为本发明中索引数组和序号数组的存储示意图。索引数组按局部排序序号(即坐标值)的升序存储各个点集中点的局部编号。序号数组按照局部编号的升序存储各个点集中点的局部排序序号。其中，局部编号指点在子集中的编号，局部排序序号指点在子集中的排序序号。

对每个子集中点的排序的流程如下：

(1)初始化j＝1；

(2)确定子集j中的第一个和最后一个点在Ndsort中的位置：

jsta＝Npoint(j)，jend＝Npoint(j)+Nstrip(j)-1；

(3)形成子集j的坐标数组(lx，ly，lz)：

lx(1～Nstrip(j))＝x(Ndsort(dsta～jend))，ly(1～Nstr(j))＝y(Ndsort(jsta～jend))，

lz(1～Nstrip(j))＝z(Ndsort(jsta～jend))；

其中，上面三个式子中的波浪线“～’表示在取值范围内逐一取值，例如1～Nstrip(j)表示在从1到Nstrip(j）的范围内逐一取值；

(4)对lx、ly和lz分别排序，形成局部索引数组（I′^x、I′^y、I′^z）和局部序号数组（R′^x、R′^y、R′^z），其中I′^x(i)存储局部排序序号为i的点的局部编号，R′^x(i)存储局部编号为i的点的局部排序序号，I′^y、I′^z和R′^y、R′^z的含义分别与I′^x和R′^x类似。

(5)分别将局部索引数组和局部序号数组转化为全局索引数组和序号数组：

I^x(jsta～jend)=I′^x(1～Nstrip(j))，I^y(jsta～jend)=I′^y(1～Nstrip(j))，

I^z(jsta～jend)=I′^z(1～Nstrip(j))，R^x(jsta～jend)=R′^x(1～Nstrip(j))，

R^y(jsta～jend)=R′^y(1～Nstrip(j))，R^z(jsta～jend)=R′^z(1～Nstrip(j))。

(6)判断j是否等于tns，如果不是则j=j+1，返回步骤（2）继续循环排序，如果是则结束排序。

步骤3、搜索给定盒子内的点

在三维空间中，以对象粒子a为中心构造边长为4h_a（h_a是现有技术中光滑粒子法中所确定的粒子a的光滑长度）的正方体盒子k。搜索给定盒子k内的点的流程如下（为叙述方便，在此仍假设条形方向为z方向）：

(1)计算包含盒子k的x和y方向边界的条形区域在x和y方向上的序号：

式中：x_lk和x_tk为盒子k的两个x方向边界的位置,y_lk和y_tk为两个y方向边界的位置。

(2)x和y方向的序号分别在Istrip_min～Istrip_max和Jstrip_min～Jstrip_max之间的条形区域对应的子集可能包含位于盒子k内的点，对这些子集按照以下程序循环搜索：

(2.1)初始化Istrip=Istrip_min，开始循环1；

(2.2)初始化Jstrip=Jstrip_min，开始循环2；

(2.3)根据公式（2）计算x和y方向的序号分别为Istrip和Jstrip的条形区域对应的子集的编号ns；

(2.4)确定子集ns中第一个和最后一个点的位置：

ista=Npoint(ns)，iend=Npoint(ns)+Nstrip(ns)-1；

(2.5)形成局部索引数组和排厅数组：

I＇^x(1～Nstrip(ns))＝I^x(ista～iend)，I＇^y(1～Nstrip(ns))＝I^y(ista～iend)，

I＇^z(1～Nstrip(ns))＝I^z(ista-iend)，R＇^x(1～Nstrip(ns))＝R^x(ista～iend)，

R＇^y(1～Nstrip(ns))＝R^y(ista～iend)，R＇^z(1～Nstrip(ns))＝R^z(ista～iend)。

(2.6)构造点列表：对局部索引数组采用现有技术中公知的二分法进行搜索，得到子集ns中在各个方向上的第一个和最后一个位于盒子内的点的局部排序序号。局部排序序号在之间的点构成x方向的点列表。类似地，根据可以分别确定y和z方向的点列表。

(2.7)寻找点列表的交集：选择点数最少的点列表，对其中的点逐个检查，判断其在其他方向的局部排序序号是否位于相应方向的第一个和最后一个点的排序序号之间。如：假设z方向的点列表中的点数最少，那么逐个检查局部编号为的点；如果且则局部编号为i＇的点在盒子k内，否则不在其中。得到点的局部编号i＇后，其全局编号可以由数组Ndst得到，即i＝Ndsort(i＇+Npoint(ns)-1)。

(2.8)判断Jstrip是否等于Jstrip_max,如果不是则Jstrip＝Jstrip+1，返回步骤(2.3)继续循环2，如果是则结束循环2，转入下一步；

(2.9)判断Istrip是否等于Istrip_max，如果不是则Istrip＝Istrip+1，返回步骤(2.2)继续循环1，如果是则结束循环1从而结束搜索。

(3)合并第(2)步中在每一个子集中搜索到的点，即为点集中位于给定盒子k内的所有点。

步骤4、条形化PIB搜索法的相邻粒子搜索

之前介绍了条形化PIB搜索法搜索给定盒子内的点的基本过程，本步骤要利用其结果进行相邻粒子的搜索，其只需要在前面过程的基础上增加少量的判断，用来判断盒子内的粒子是否与对象粒子a构成相邻粒子对。

由于现有技术中光滑粒子法中的光滑函数具有紧支性，在对某个对象粒子进行近似计算时，只需要位于其支持域内的粒子，这些粒子称为该粒子的相邻粒子。为了示意方便，图3示出了在二维空间中对象粒子的相邻粒子示意图，其中对象粒子a的支持域是以其为中心、以2h_a为半径的圆，类似地，在三维空间中，对象粒子a的支持域是以其为中心、以2h_a为半径的球。

从严格意义上讲，支持域中任意粒子b构成对象粒子a的相邻粒子的条件是（此处a和b均为粒子编号）：

d_ab<2h_a    （4）

以上公式（4）中d_ab指粒子a、b间的距离。按公式（4）定义的相邻粒子是不对称的，即粒子b为对象粒子a的相邻粒子时，对象粒子a不一定为粒子b的相邻粒子，反之亦然。因此，在实际计算中，经常基于对称光滑长度来定义相邻粒子，这样使光滑粒子法的方程完全对称，从而获得更精确、稳定的结果。采用对称光滑长度定义时，成为相邻粒子需满足的条件是：

$d_{\mathrm{ab}}<2{\stackrel{\&OverBar;}{h}}_{\mathrm{ab}}=(h_a+h_b)---\left(5\right)$

如果粒子a和b满足上述公式（5）的条件，则它们互为彼此的相邻粒子，构成一个相邻粒子对。在本发明中即采用对称光滑长度定义相邻粒子。

判断盒子内的粒子与对象粒子a是否作为相邻粒子对的方法为：如果搜索到的粒子b满足：

a<b，d_ab<2h_a， $d_{\mathrm{ab}}<2{\stackrel{\&OverBar;}{h}}_{\mathrm{ab}}---\left(6\right)$

或者

a>b，d_ab≥2h_b， $d_{\mathrm{ab}}<2{\stackrel{\&OverBar;}{h}}_{\mathrm{ab}}---\left(7\right)$

则粒子a和粒子b记录为相邻粒子对。式（6）和（7）在式（5）的相邻粒子对定义的基础上附加了条件，其目的是为了防止重复记录相邻粒子对而引起错误，即：在搜索粒子a的相邻粒子时将粒子a和粒子b记录为相邻粒子对，而在搜索粒子b的相邻粒子时又再次将他们作为相邻粒子对。

下面通过两个具体实施例来测试条形化PIB搜索法的性能，并将其与PIB搜索法、树形搜索法进行比较。在具体实施例中，所有计算时间的统计结果均得自同一计算机，所用的计算机的主要配置为：1.40GHz Intel CPU，1.89GB内存。

实施例1、条形化PIB搜索法的性能测试及其与PIB搜索法比较。

图4、图5为二维测试使用的两种点集：一种为条形分布，y方向的点数固定为n_y=5，x方向的点数n_x可变且远大于n_y；另一种为正方形分布，x方向和y方向的点数均为n。对这两个点集中的每个点，定义以点为中心、尺寸为Δ_b＝6Δ（Δ为粒子间距）的方形盒子。搜索每个盒子内的点，以此来测试搜索算法的效率。

首先，研究条形尺寸对条形化PIB搜索法的效率的影响。选用方形分布的点集，取点数为250000（n＝500）。采用13种不同条形的条形尺寸Δ_s进行搜索计算，Δ_s/Δ_b由0.125变化到100。下表给出了在不同Δ_s/Δ_b值下的搜索时间。由该表可见，搜索时间随着Δ_s/Δ_b的增加先减少，到0.5时达到最少，之后开始单调增加。当Δ_s/Δ_b=100时，只有一个子集，条形化PIB搜索法退化为PIB搜索法，搜索时间为12.05s。从下表还可以看到，Δ_s/Δ_b在0.125到10（跨越将近两个数量级）间变化时，搜索时间仅在1.21～2.25s间变化，与PIB搜索法的相对时间在10%～19%间变化。这说明，条形化PIB搜索算法在相当宽范围的条形尺寸下都可以获得很高的效率，所以它能够适应盒子尺寸具有大的不一致性的问题。根据上述结果，本发明采用以下方法确定条形尺寸Δ_s：

${\mathrm{\&Delta;}}_s=\left\{\begin{array}{c}0.5{\mathrm{\&Delta;}}_b\\ 0.5\sqrt{{\mathrm{\&Delta;}}_{b\max }{\mathrm{\&Delta;}}_{b\min }}\end{array}\right.---\left(8\right)$

上式中，上面一行结果适用于Δ_b一致的问题，而下面一行结果适用于Δ_b不一致的问题，其中Δ_bmax和Δ_bmin分别为最大和最小的盒子尺寸。

然后，研究不同点数和点分布情况对条形化PIB搜索法的效率的影响，并将其与PIB搜索法进行比较。图6为针对条形分布和方形分布的点集在不同点数情况下的搜索耗时。由该图可见，条形化PIB搜索法对于两种点集的搜索效率都很高，差别不大。而PIB搜索法对于两种点集的搜索效率差别很大，对于条形分布的点集，其效率略低于条形化PIB搜索法；对于方形分布的点集，随着点数的增多，其效率急剧降低，远低于条形化PIB搜索法。综上可见，条形化PIB搜索法的性能明显地优于PIB搜索法。

图7、图8为三维测试使用的点集。与二维测试类似，仍考虑两种不同分布的点集：一种为条形分布的点集,x和y方向的点数固定为n_x=n_y=5，z向的点数n_z远大于n_x和n_y；另一种为立方形分布，各个坐标方向的点数均为n。三维测试中的盒子为以每个点为中心、边长为6倍点间距的立方形盒子。

图9比较了在不同点数和点分布情况下条形化PIB搜索法与PIB搜索法的搜索时间。由该图可见，对于条形分布的点集，条形化PIB搜索法与PIB搜索法的搜索时间很接近；对于立方形分布的点集，前者远远小于后者。另外，由该图可见，PIB搜索法的搜索时间受点集的分布影响很大，而条形化PIB搜索法的搜索时间受点集的分布影响很小。图10、图11给出在不同点数情况下PIB搜索法与条形化PIB搜索法的搜索耗时比。由该图可见，对于条形分布的点集，搜索耗时比约为1，且基本不随点数的增加而变化，这说明对于条形分布的点集，条形化PIB搜索法的效率与PIB搜索法相当。对于方形分布的点集，搜索耗时比随点数的增加而迅速增大。当点数为300000时，搜索耗时比达到25.6，这说明对于立方形分布的点集，条形化PIB搜索法远远优于PIB搜索法。

实施例2、条形化PIB搜索法与树形搜索法比较

首先，采用图7中的均匀立方形分布的粒子集来比较条形化PIB搜索法与树形搜索法的相邻粒子搜索效率。粒子的光滑长度设置为1.5倍粒子间距，每个粒子的相邻粒子数平均约为110个。图11比较了在不同粒子数的情况下条形化PIB搜索法和树形搜索法完成相邻粒子搜索的耗时。由该图可见，条形化PIB搜索算法所需的时间小于树形搜索算法，这说明前者的效率高于后者。图12为树形搜索法与条形化PIB搜索法的搜索耗时比。由该图可见，搜索耗时比在3.6～4.6的范围内变化，这说明条形化PIB搜索法的搜索速度为树形搜索法的4倍左右。

然后，考虑非均匀粒子分布情况。以模拟卵形钢弹以30°倾角斜侵彻铝板的三维侵彻模型作为测试问题。图13为计算采用的三维光滑粒子法模型。为节约计算时间，沿弹体和靶体的对称面分割，采用半模型进行计算。在对称面处应用镜像虚粒子来施加约束。弹体的初始粒子间距为1mm；靶体在中心区域的初始粒子间距为Δ=1mm，在外围区域为Δ=2mm；模型共计使用202308个粒子（不含虚粒子），其中弹体粒子数为6892，靶体粒子数为195416。计算中，粒子的初始光滑长度设为h＝1.5Δ。图14为采用条形化PIB搜索法和树形搜索法进行相邻粒子搜索时一个三维光滑粒子法的计算循环的耗时。由该图可见，采用树形搜索法进行相邻粒子搜索时，一个计算循环的耗时为71.0s，其中相邻粒子搜索时间为48.3s，占总耗时的68%；而采用条形化PIB搜索法进行相邻粒子搜索时，一个计算循环的耗时减小为34.4s，其中相邻粒子搜索时间仅为12.7s，占总时间的37%。由此可见，在非均匀粒子分布的情况下，条形化PIB搜索法的相邻粒子搜索效率仍为树形搜索法的4倍左右。相比于树形搜索法，它能使光滑粒子法的整体效率提高约一倍。类似地，本发明中的条形化PIB搜索法也可以用于汽车碰撞、飞机降落以及爆炸冲击等工程实际问题的模拟仿真中。

根据上述具体实施例的介绍可知，本发明具有如下显著的技术效果：

（1）本发明提出的条形化PIB搜索方法是在PIB搜索方法的基础发展起来的，克服了PIB搜索法对粒子点分布的依赖，具有比PIB搜索法和树形搜索法更高的计算效率；

（2）条形化PIB搜索方法适用于各种商用软件，如：ANSYS、LS-DYNA、ABAQUS等，可嵌入软件中提高相邻粒子、接触搜索算法的计算效率，提高软件模拟大规模工程问题的能力。

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

Claims (9)

1.一种模拟仿真技术中的相邻粒子搜索方法，其特征在于，所述方法包括以下步骤： 

步骤1：将点集划分为一系列条形的子集； 

步骤2：分别对每个子集内的点排序； 

步骤3：搜索给定盒子内的点； 

步骤4：判断盒子内的粒子是否与对象粒子形成相邻粒子对，从而搜索到所有相邻粒子； 

其中，在所述步骤1中，将点集占据的最小区域在垂直于条形方向的方向上进行分割，以得到若干小的条形区域，每一个小的条形区域中的点即构成为一个子集，其中条形方向是点集宽度最大的方向，即坐标跨度最大的方向。 

2.根据权利要求1所述的方法，其特征在于，在所述步骤2中，对每个子集中的点按各个坐标方向的坐标值排序，点排序的结果是索引数组和序号数组，其中该索引数组按局部排序序号的升序存储各个点集中点的局部编号，该序号数组按照局部编号的升序存储各个点集中点的局部排序序号，其中局部编号指点在子集中的编号，局部排序序号指点在子集中的排序序号。 

3.根据权利要求1所述的方法，其特征在于，在所述步骤3中，先确定可能包含位于盒子内的点的子集，再利用PIB搜索法逐个确定这些子集中位于盒子内的点，最后将在每个子集中搜索到的点合并即得到盒子内的所有点。 

4.根据权利要求1所述的方法，其特征在于，在所述步骤1中，当三维空间的条形区域的尺寸为Δ_s、条形的方向为z方向时，将点集占据的区域在x和y方向上进行分割，其分割数分别为： 

tns^x＝int[(x_max-x_min)Δ_s]+1和tns^y＝int[(y_max-y_min)Δ_s]+1   (1) 

公式(1)中：x_max、x_min分别为点的最大和最小x坐标，y_max、y_min分别为点的最大和最小y坐标； 

将点集占据的区域在x和y方向上分割后，得到tns＝tns^x×tns^y个小的条形区域，对每个条形区域按如下方式编号： 

ns(ns^x,ns^y)＝(ns^y-1)tns^x+ns^x   (2) 

公式(2)中：ns为条形区域的编号，ns^x和ns^y分别为条形区域在x和y方向的序号； 

对于任意给定的点i，该点所在的条形区域在x和y方向上的序号按下式计算： 

ns_i ^x＝min(int[(x_i-x_min)Δ^x _s]+1,tns^x) 

                                         (3) 

ns_i ^y＝min(int[(y_i-y_min)Δ^y _s]+1,tns^y) 

公式(3)中：和分别为x和y方向的精确的条形尺寸； 

将按公式(3)计算的序号代入公式(2)，即得到点i的所在的条形区域的编号ns_i，采用上述方法确定每一点所在的条形区域的编号，所有在同一个条形区域中的点构成一个子集，子集的编号为所在条形区域的编号。 

5.根据权利要求4所述的方法，其特征在于，在所述步骤1中，在划分子集时构造以下数组来记录子集的相关信息从而方便后续的搜索操作：长度为总点数N的数组Strip和Ndsort，其中Strip存储每一个点的所在的子集编号，Ndsort按点所在子集的编号的升序存储点的编号；长度为tns的数组Nstrip和Npoint，其中Nstrip记录每一个子集中的点数，Npoint记录每一个子集中的第一个点在Ndsort中的位置，其中，构造这些数组的流程如下： 

(1)初始化数组：Nstrip＝0； 

(2)根据公式(2)确定每一个点i所在的子集编号ns_i； 

(3)存储点i的子集编号到数组Strip：Strip(i)＝ns_i； 

(4)修改子集ns_i中的点数：Nstrip(ns_i)＝Nstrip(ns_i)+1； 

(5)计算每一个子集j中的第一个点在Ndsort中的位置： 

Npoint(1)＝1,Npoint(j)＝Npoint(j-1)+Nstrip(j-1)； 

(6)初始化数组：Nstrip＝0； 

(7)将点按其所在的子集的编号的升序存储到Ndsort中： 

Ndsort(Nstrip(Strip(i))+Npoint(Strip(i)))＝i，Nstrip(Strip(i))＝Nstrip(Strip(i))+1。 

6.根据权利要求5所述的方法，其特征在于，在所述步骤2中，对每个子集中点的排序的具体流程如下： 

(1)初始化j＝1； 

(2)确定子集j中的第一个和最后一个点在Ndsort中的位置： 

jsta＝Npoint(j),jend＝Npoint(j)+Nstrip(j)-1； 

(3)形成子集j的坐标数组(lx，ly，lz)： 

lx(1～Nstrip(j))＝x(Ndsort(jsta～jend))，ly(1～Nstrip(j))＝y(Ndsort(jsta～jend))， 

lz(1～Nstrip(j))＝z(Ndsort(jsta～jend))； 

其中，上面三个式子中的波浪线“～”表示在取值范围内逐一取值； 

(4)对lx、ly和lz分别排序，形成局部索引数组(I^'x、I^'y、I^'z)和局部序号数组(R^'x、R^'y、R^'z)，其中I^'x(i)存储局部排序序号为i的点的局部编号，R^'x(i)存储局部编号为i的点的局部排序序号，I^'y、I^'z和R^'y、R^'z的含义分别与I^'x和R^'x类似； 

(5)分别将局部索引数组和局部序号数组转化为索引数组和序号数组： 

I^x(jsta～jend)＝I^'x(1～Nstrip(j))，I^y(jsta～jend)＝I^'y(1～Nstrip(j))， 

I^z(jsta～jend)＝I^'z(1～Nstrip(j))，R^x(jsta～jend)＝R^'x(1～Nstrip(j))， 

R^y(jsta～jend)＝R^'y(1～Nstrip(j))，R^z(jsta～jend)＝R^'z(1～Nstrip(j))；

(6)判断j是否等于tns，如果不是则j＝j+1，返回步骤(2)继续循环排序，如果是则结束排序。 

7.根据权利要求6所述的方法，其特征在于，在所述步骤3中，以对象粒子a为中心构造边长为4h_a的正方体盒子k，其中h_a是粒子a的光滑长度， 

搜索给定盒子k内的点的具体流程如下： 

(1)计算包含盒子k的x和y方向边界的条形区域在x和y方向上的序号： 

上式中：x_lk和x_tk为盒子k的两个x方向边界的位置,y_lk和y_tk为两个y方向边界的位置； 

(2)x和y方向的序号分别在Istrip_min～Istrip_max和Jstrip_min～Jstrip_max之间的条形区域对应的子集可能包含位于盒子k内的点，对这些子集中的点进行搜索； 

(3)合并上述步骤(2)中在每一个子集中搜索到的点，即得到点集中位于给定盒子k内的所有点。 

8.根据权利要求7所述的方法，其特征在于，在所述步骤3的上述步骤(2)中，对这些子集按照以下方式搜索： 

(2.1)初始化Istrip＝Istrip_min，开始循环1； 

(2.2)初始化Jstrip＝Jstrip_min，开始循环2； 

(2.3)根据公式(2)计算x和y方向的序号分别为Istrip和Jstrip的条形区域对应的子集的编号ns； 

(2.4)确定子集ns中第一个和最后一个点的位置： 

ista＝Npoint(ns)，iend＝Npoint(ns)+Nstrip(ns)-1； 

(2.5)形成局部索引数组和排序数组： 

I^'x(1～Nstrip(ns))＝I^x(ista～iend)，I^'y(1～Nstrip(ns))＝I^y(ista～iend)， 

I^'z(1～Nstrip(ns))＝I^z(ista～iend)，R^'x(1～Nstrip(ns))＝R^x(ista～iend)， 

R^'y(1～Nstrip(ns))＝R^y(ista～iend)，R^'z(1～Nstrip(ns))＝R^z(ista～iend)； 

(2.6)构造点列表：对局部索引数组采用二分法进行搜索，得到子集ns中在各个方向上的第一个和最后一个位于盒子内的点的局部排序序号(和)；局部排序序号在之间的点构成x方向的点列表，同样根据和分别确定y和z方向的点列表； 

(2.7)寻找点列表的交集：选择点数最少的点列表，对其中的点逐个检查，判断其在其他方向的局部排序序号是否位于相应方向的第一个和最后一个点的排序序号之间：当z方向的点列表中的点数最少，则逐个检查局部编号为的点；如果 且则局部编号为i'的点在盒子k内，否则不在其中；得到点的局部编号i'后，其全局编号由数组Ndsort得到，即i＝Ndsort(i'+Npoint(ns)-1)； 

(2.8)判断Jstrip是否等于Jstrip_max，如果不是则Jstrip＝Jstrip+1，返回步骤(2.3)继续循环2，如果是则结束循环2，转入步骤(2.9)； 

(2.9)判断Istrip是否等于Istrip_max，如果不是则Istrip＝Istrip+1，返回步骤(2.2)继续循环1，如果是则结束循环1从而结束步骤(2)的搜索。 

9.根据权利要求8所述的方法，其特征在于，在所述步骤4中，以如下方式来判断盒子内的粒子是否与对象粒子a形成相邻粒子对： 

如果搜索到的粒子b满足： 

或者 

则粒子a和粒子b记录为相邻粒子对，其中a和b均为粒子编号，上式中d_ab指粒子a、b间的距离，其中h_a是粒子a的光滑长度，h_b是粒子b的光滑长度。