CN106485030B - 一种用于sph算法的对称边界处理方法 - Google Patents

Abstract

本发明公开了一种用于SPH算法的对称边界处理方法，其特征在于，包括如下步骤：创建待进行仿真分析的仿真模型并进行初始化设置；计算仿真模型所对应的每一粒子的压力值、应力张量以及人工粘性值；在所设定的每一个时间步内对所述仿真模型进行搜索分析；计算仿真模型所对应的状态方程以及控制方程，若仿真模型存在对称边界，则对仿真模型进行作用对信息增补处理；基于所设定终止条件，输出对应的仿真结果。本发明通过对现有作用对信息进行筛查增补，使得本发明不需要额外生成镜像粒子、也不需要再对镜像粒子划分背景网格进行搜索的同时，去掉了额外生成镜像粒子和搜索的时间消耗，实现了将对称边界处理和正常的搜索计算统一处理。

Description

一种用于SPH算法的对称边界处理方法

技术领域

本发明涉及计算机仿真技术领域，具体的说是涉及一种用于SPH算法的对称边界处理方法。

背景技术

通常我们关心的物体运动和变形的规律可以采用微分方程来描述，称这种微分方程叫做物理控制方程，实际上大部分的运动物体都是形状复杂的，物理控制方程必须借助于计算机帮助进行数值计算，也叫数值仿真或数值模拟，其意思是根据客观规律采用数值计算的方法模拟物体在一定条件下的运动和变形情况，这属于一种预测，现在技术条件下的预测结果与实际情况已经能够非常接近。因而，数值仿真计算在科学研究和工程应用中已经成为并列于实验手段和理论手段的第三种研究方法，其具有结果直观，可重复性高，条件要求低和成本低等优势。

通常，数值计算方法可以分为两种：一种是基于网格的数值计算方法，另一种则是无网格数值计算方法；其中基于网格的数值计算方法分为拉格朗日方法以及欧拉方法：所述拉格朗日方法是让网格随着物体的运动和变形而变化，该方法主要用在固体的数值计算，其能够精确捕获固体的形状变化，对于不规则外形可以用不规则的网格分布来描述，计算效率高，适应性强，广泛用于固体力学计算中，在机械行业中也有大量的应用；但是正因为拉格朗日方法的网格随着变形而变化，当变形很大的时候，网格变形很严重，此时就不能满足有网格数值计算方法的基本要求，因此会造成时间步长减少，每一次计算的耗时增加，甚至因为网格畸变严重导致计算失败，限制了这种方法应用；所述欧拉方法是让网格固定不变，每一次获得的结果都是空间固定位置的物理量，主要用于流体问题的计算，因为网格固定不变，因此不会因为网格变形影响计算，但是正是因为固定的网格，因此其不能描述固体的变形，并且在用于流体计算中没有确认多种流体之间界面的措施，因此需要额外提出确认流体界面的方法，欧拉方法目前尚没有像拉格朗日方法那么有效的自然处理界面的方法，因此其在流体中采用的处理界面的方法只能在一定程度上满足计算的要求，并不十分精确，基本能够满足流体计算要求，但是对于固体计算来说让不能满足固体边界的精确描述。

从两种有网格方法的特点来看，有网格方法都存在一些不适用的缺陷，虽然可以相互补充，但由于描述方法的不同而不能一起使用，如果要满足固体力学中的数值计算要求，就需要去除网格，满足流体力学计算的要求，就需要一种拉格朗日方法描述多物质流体的界面。

而基于无网格方法的SPH方法刚好可以满足这些要求，SPH方法全称为光滑粒子流体动力学方法，是一种拉格朗日格式的无网格方法，其既可以用于固体力学数值计算，也可以用于流体力学数值计算；其在用于固体力学时，因为没有网格，因此网格畸变的问题不存在于SPH方法，其用于流体力学时，因为是拉格朗日方法，因此不需要额外的界面追踪技术就能精确描述多种流体的交界面。

随着研究的深入，SPH方法被逐步应用于水下爆炸和高速碰撞，在海洋波浪模拟和液面晃动的模拟中也有成功应用，近年来被广泛用于航天器空间碎片防护的研究。但是随着SPH方法逐步被用于实际的工程应用中，现有的SPH方法也存在一定的缺陷，其具体包括下述问题：

在现有的SPH方法中尚没有出现在对固壁边界进行处理时，特别是对称边界进行处理时，能够高效进行处理的方法。关于固壁边界的处理主要有两种，第一种是在边界上预先布置一层粒子，所有粒子与这层例子都有一个反作用力，当粒子靠近这层粒子的时候就会受到相反方向的作用，组织粒子穿透，当粒子远离的时候就会受到拉伸作用，防止分离，但是这种方法弊端过多，实用性不高，所以出现了第二种方法-镜像粒子方法，镜像粒子方法就是在对称的边界的另一边布置与这边呈镜面对称的粒子，因为如果模型是对称的，那么完整模型计算的结果也是这样的，因此只要每一步做镜像处理，就等于只计算一半的粒子就可以获得所有粒子的状态，这种方法的好处是非常精确，与完整的计算一样，但坏处是需要额外生成镜像粒子，并且镜像粒子数量与实际粒子的数量相同，并且需要再做一次搜索，造成资源的不必要浪费。

发明内容

鉴于已有技术存在的缺陷，本发明的目的是要提供一种用于SPH算法的对称边界处理方法，其不需要额外生成镜像粒子、也不需要再对镜像粒子划分背景网格进行搜索的同时，去掉了额外生成镜像粒子和搜索的时间消耗，实现了将对称边界处理和正常的搜索计算统一处理。

为了实现上述目的，本发明的技术方案：

一种用于SPH算法的对称边界处理方法，其特征在于，包括如下步骤：

步骤1、创建待进行仿真分析的仿真模型并对所创建的仿真模型进行初始化设置；

步骤2、计算经初始化设置后的仿真模型所对应的每一粒子的压力值、应力张量以及人工粘性值；

步骤3、在所设定的每一个时间步内对所述仿真模型进行搜索分析，以确定每一粒子所对应的近邻粒子；

步骤4、计算仿真模型所对应的状态方程以及控制方程，在进行计算时，若仿真模型存在对称边界，则对仿真模型进行作用对信息增补处理；

步骤5、基于所设定终止条件，输出对应的仿真结果。

进一步的，作为本发明的优选

若所述对仿真模型为1/2对称模型，所述对仿真模型进行作用对信息增补处理包括：

步骤41、分别确认每一分区是否存在对称边界，若当前的分区存在对称边界，则对当前的分区中所对应的所有作用对进行循环，并判断某一作用对中是否存在至少一个粒子到对称边界的距离小于自身紧支域，是则确认可能存在具有镜像作用的粒子对；

步骤42、进一步确定该作用对的镜像关系，即判断该作用对中的任意一种粒子的镜像粒子是否是另一粒子的作用对，所述的判断策略为：获取该作用对中的任意一种粒子的镜像粒子的镜像位置，并判断所获取的镜像位置与另一粒子的间距与两个粒子紧支域之和的差值的绝对值是否小于所确定的判断阈值，是则确定该作用对为具有镜像作用的粒子对即确认该作用对中的任意一个粒子的镜像粒子与另一粒子互为作用对，且若步骤41所述的某一作用对中存在其到分区边界的距离小于自身紧支域，则确认该粒子的镜像粒子与该粒子互为作用对；依次对所有作用对循环，以获得全部的关于对称边界的作用对信息，对仿真模型进行作用对信息增补处理。

进一步的，作为本发明的优选

若所述对仿真模型为多重对称模型即1/4对称模型或者1/8对称模型，所述对仿真模型进行作用对信息增补处理包括：

步骤41′、分别独立确认每一并行分区是否存在对称边界，若当前的并行分区存在对称边界且对称边界个数至少为2个，则分别沿着x方向、y方向以及z方向对当前的并行分区中所对应的所有作用对进行循环，并在依次对上述3个方向进行循环的过程中，判断当前的作用对中是否存在至少一个粒子到其中任意一个对称边界的距离均小于自身紧支域，是则确认可能存在具有镜像作用的粒子对；

步骤42′、进一步确定该作用对的镜像关系，因为是对称的关系，所以需要考虑一个粒子在各个对称面上的镜像是否均是另一个粒子的作用对，即判断该作用对中的任意一种粒子的在各个对称面上的镜像粒子是否均是另一粒子的作用对，所述的判断策略为：首先获取该作用对中的任意一种粒子在某一对称面上镜像粒子的镜像位置，并判断所获取的镜像位置与另一粒子的间距与两个粒子紧支域之和的差值的绝对值是否小于所确定的判断阈值，是则确定确认该粒子在当前的对称面上镜像粒子与另一粒子互为作用对，且若该作用对中存在其到分区边界的距离小于自身紧支域，则确认该粒子在当前的对称面上的镜像粒子与该粒子互为作用对；其次获取该作用对中的任意一种粒子在其他对称面上镜像粒子的镜像位置，并判断所获取的镜像位置与另一粒子的间距与两个粒子紧支域之和的差值的绝对值是否小于所确定的判断阈值，是则确定确认该粒子在所述对称面上镜像粒子与另一粒子互为作用对，且若步骤41′所述当前的作用对中存在其到分区边界的距离小于自身紧支域，则确认该粒子在所述对称面上的镜像粒子与该粒子也互为作用对；依次对所有作用对循环，以获得全部的关于对称边界的作用对信息，对仿真模型进行作用对信息的进行增补；同时所述作用对信息包括相互作用的粒子的标号以及各自所对应的作用方向，并在获取作用对信息之后，根据对称记号计算人工粘性和相互作用函数，同样记录在作用对信息之中。

与现有技术相比，本发明的有益效果：

在对对称边界进行处理时，本发明通过对现有作用对信息进行筛查增补，使得本发明不需要额外生成镜像粒子、也不需要再对镜像粒子划分背景网格进行搜索的同时，去掉了额外生成镜像粒子和搜索的时间消耗，实现了将对称边界处理和正常的搜索计算统一处理。

附图说明

图1为本发明所述的用于SPH算法的对称边界处理方法对应的流程步骤图；

图2为本发明所述的用于SPH算法的对称边界处理方法对应的流程步骤例图；

图3为本发明所述计算域被划分为均等的4块区域的例图；

图4为本发明所述搜索域的大小大于当前待搜索的区域大小的例图；

图5a为对图4进行有向搜索时所对应的例图1；

图5b为对图4进行有向搜索时所对应的例图2；

图5c为对图4进行有向搜索时所对应的例图3；

图5d为对图4进行有向搜索时所对应的例图4。

具体实施方式

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

如图1-图2所示，本发明所述一种用于SPH算法的对称边界处理方法，其特征在于，包括如下步骤：

(1)、创建待进行仿真分析的仿真模型并对所创建的仿真模型进行初始化设置；由于SPH仿真计算过程的前处理可以与有限元仿真计算兼容，即可以通过网格来确定其初始的粒子分布情况，即需要完成粒子的分布设定，模拟条件的设定，粒子所代表材料的属性和粒子的初始状态等设置。具体的，作为本发明的优选实例，所述的步骤1其包括：

步骤101、创建仿真模型并设定所创建的仿真模型所对应的网格文件，获得对应的网格信息并存入临时的网格数据中，所述网格信息包括网格序号以及网格材料号，且网格信息包含与其所对应的节点号；

步骤102、设定所创建的仿真模型所对应的节点文件，获得对应的节点信息并存入临时的节点数据中，所述节点信息包括节点号以及节点坐标；

步骤103、结合网格信息与节点信息计算仿真模型所对应的各个网格的中心坐标以及每一网格所对应的体积数据，并基于所计算的体积数据计算各自对应的等效直径；

步骤104、将上述各网格转换为所对应的粒子，其中每一粒子的位置即其所对应的网格的中心坐标，其大小即其所对应的网格的等效直径，其材料号即其所对应的网格材料号；

步骤105、设定材料文件，获取仿真模型所对应的每种材料的密度，泊松比，弹性模量以及材料模型的类型和参数等数据；

步骤106、根据每一粒子的所对应的网格材料号将所设定的材料的属性和模型参数赋给每一粒子；

步骤107、设定边界条件和初始条件，获取仿真模型所对应的条件值和条件作用的组件号，存在临时的初边值数据中；

步骤108、设定组件信息文件并获得仿真模型所对应的组件数，若按顺序打开组件文件，则依次设定对应的组件号、组件类型；进一步的，若当前的组件号和组件类型满足当前所对应的初边界条件，则继续设定当前组件的成员，获取相关网格的序号，将初边值条件赋给网格相应的粒子；若边界条件作用于节点，则获取该节点属于的单元数据，将边界条件赋给对应单元转换得到的粒子；

步骤109、设定初始控制参数；所述计算控制参数包括当前时间、总体计算时长、当前时间步长、总体循环数、当前时间步数、最小时间步长、最大时间步长、对称信息和保存结果设定、并行信息等信息。

步骤1010、设定每种材料的属性和模型参数，将所设定的材料参数和模型参数按照粒子的材料号赋给每个粒子。

(2)、计算经初始化设置后的仿真模型所对应的每一粒子的压力值、应力张量以及人工粘性值；具体的，作为本发明的优选实例，所述的压力数据的计算过程为：首先对每一粒子进行循环并获得每一粒子的密度和比内能，将所获得的数据作为计算压力的输入参数；其次根据各粒子所赋予的状态方程中确定其压力计算函数的类型；最后基于所确定的全部输入参数以及压力计算函数，计算每个粒子的压力后结束粒子循环过程；具体的，作为本发明的优选实例，所述的应力值的计算过程为：首先对每一粒子进行循环并获得每一粒子的应变率、旋转率以及应力；其次基于应力率公式计算出每一粒子所对应的应力率；再次确定当前的时间步长后，用时间步长乘以所对应的应力率得到当前所计算的粒子的应力张量，并通过将该粒子的应力加上应力张量得到临时的更新后的应力，将该粒子的应力带入所设定的屈服条件，判断是否屈服，如果屈服则用径向返回方法获得屈服后的应力，如果没有屈服则不进行处理；最后将判断后的应力作为该粒子更新后的应力并结束粒子循环过程；具体的，作为本发明的优选实例，所述人工粘性值的计算过程为：首先对每个作用对进行循环并依次选定其中一个作用对中的粒子i和粒子j，计算这一对作用粒子相互之间的作用力，即粒子i和粒子j的相对坐标(x_i-x_j，y_i-y_j，z_i-z_j)和相对速度(v_xi-v_xj，v_yi-v_yj，v_zi-v_zj)，其中x_i，y_i和z_i表示粒子i的x方向，y方向和z方向坐标，x_j，x_j和z_j表示粒子j的x方向，y方向和z方向坐标，v_xi，v_yi和v_zi表示粒子i的x方向，y方向和z方向速度，v_xj，v_yj和v_zj表示粒子j的x方向，y方向和z方向速度；随后粒子i的光滑长度h_i，密度rho_i，声速c_i和j粒子的光滑长度h_j，密度rho_j，声速c_j；最后计算出这一对作用粒子相互之间的作用力，并将该作用力作为这一个作用对的人工粘性力，并结束作用对循环。

(3)、在所设定的每一个时间步内，对所述仿真模型进行搜索分析，以确定每一粒子所对应的近邻粒子；如图2(图中虚线表示不是必要步骤，其可以自由选择添加该步骤)，所述的搜索分析包括步骤31确定仿真模型所对应的各并行分区以及各并行分区边界以及步骤32对每一并行分区进行区域搜索；

具体的，作为本发明的优选实例，所述步骤31包括如下步骤：

步骤311、首先对仿真模型所对应的所有粒子进行循环，分别找到在三维空间的X方向、Y方向、Z方向上，坐标最小所对应的粒子I以及坐标最大所对应的粒子II；其次将粒子II的坐标与粒子I的坐标相减，得到X、Y、Z三个方向上的最大跨度值；最后通过将所述仿真模型所对应的每一粒子的坐标依次与粒子I的坐标相减并除以所获得的该坐标方向上的最大跨度值，确定出每一粒子的相对无量纲坐标(X_j-X_i/X₀,Y_j-Y_i/Y₀,Z_j-Z_i/Z₀),其中，X_j、Y_j、Z_j分别表示X方向、Y方向、Z方向上的粒子坐标；X_i、Y_i、Z_i分别表示X方向、Y方向、Z方向上坐标最小的粒子I的坐标，X₀表示X方向上的最大跨度值，Y_j表示Y方向上的最大跨度值，Z₀表示Z方向上的最大跨度值，j≥1；

步骤312、分别计算所有粒子所对应的无量纲坐标在X方向、Y方向、Z方向上的平方和，并将平方和最小的无量纲坐标所构成的向量方向确定为主方向；

步骤313、基于所确定的主方向，按照所有粒子与主方向的相对坐标从小至大进行排序：进一步的，作为本发明的优选例，所述步骤313采用浮点数的桶排序方法对所有粒子进行排序，下述以具体实例为例进行说明：首先获取所有粒子在主方向上的相对无量纲坐标值，将所获得相对无量纲坐标值进行坐标格式化后转换为相应的字符串，由于每个转化后的字符串具有相同的长度，因此所有相对无量纲坐标均不会为负值，且每一个字符串顺序为第1位是0，第2位是小数点，第3位至第8位是浮点数小数点后6位数值，第9位是e，即表示指数开始，第10位是指数的符号，其为负值，第11到12位是指数的十位和个位；其次，设置m个桶(本例设置0到9号，一共9个桶)后，先对各字符串的小数位进行排序：依次对每一位小数位循环时，将各粒子放入对应的序号的桶中，且按所设定的顺序将先搜索到的粒子先放进去即如果对第1位小数循环，当前粒子的第一位小数位4的放入四号桶，完成第一位小数位的循环后，按照0到9的顺序将每个桶中的粒子按序输出(就是说从0号桶开始输出，先放进0号桶的先输出)以形成一个新的粒子序列，然后以这个新序列开始对下一位小数位进行排序，直到完成所有小数位的排序；随后以前述完成小数位排序所对应的输出序列开始对指数位进行排序，但是此时的排序方法应与小数相反，这是由于指数位的符号是负号，通常越大的指数表示这个数越小，所以对每一个指数位循环将各粒子按序放入桶时的对应关系为：首先若当前排序位的值为n，则将粒子放入9-n号桶中；然后按照放入的顺序将各桶中的粒子输出，此时所获得的序列作为下一个指数位排序的初始序列，直至完成所有指数位排序，此时所获得的序列就是全部粒子按照主方向相对坐标从小到大的粒子序列；

步骤314、确定待划定的并行分区数量并确定每一并行分区所对应的初始粒子数以及各并行分区边界位置所对应的分界粒子序号；进一步的，作为本发明的优选所述步骤314包括：步骤3141、获取确定仿真模型所对应的并行分区数量；如依据当前仿真环境即当前仿真计算机(特别是cpu)的配置参数确定并行分区数量；步骤3142、按照仿真模型所对应的所有粒子数与所确定的并行分区数量相除所获得的整数加一设定每一并行分区所对应的初始的粒子数；步骤3143、确定各并行分区边界位置所对应的分界粒子序号，即按照第一个并行分区与第二并行分区边界位置所对应的分界粒子序号为第一个并行分区所具有的粒子数，且第n-1个并行区与第n个并行区的分界粒子序号为第n-2个并行分区与第n-1个并行分区的分界粒子序号加上第n-1个并行区的粒子数的策略，对所有的并行分区进行循环，以确定各并行分区边界位置所对应的分界粒子序号，其中n≥3，如第二个并行分区与第三个并行分区的分界粒子序号为第一个并行分区与第二个并行分区的分界粒子序号加第二个并行区所具有的粒子数；

步骤315、确定每一并行分区边界位置；具体的，作为本发明的优选实例，所述步骤315包括：按照步骤313排好的序列进行遍历，当遍历到某一个并行分区的分界粒子序号时，提取该分界粒子序号所对应的粒子的编号，并找到该编号对应的粒子在主方向上的绝对坐标，此绝对坐标即为当前并行分区与下一个并行分区的分界位置；完成所有粒子序列的遍历以得到所有并行分区的上边界和下边界，其中上边界即为当前的并行分区与前一并行分区的边界坐标，下边界即为其与后一并行分区的边界坐标；且第一个并行分区上边界是主方向起始位置所对应的粒子坐标，最后一个并行分区下边界是主方向终止位置所对应的粒子坐标；

由于各并行分区之间的相对关系是明确的、有序的，因此计算开始的时候通过上述方法确定各并行分区后，在后续计算中只需要对每一个并行分区内部的粒子进行排序，然后通过各并行分区通信来确定各个并行分区区的粒子数量变化，根据数量变化调整分界位置，即所述步骤316所述的内容；

步骤316、划定并行分区并对每一并行分区的粒子进行排序，并通过各并行分区通信来确定各个并行分区内的粒子数量变化情况，并基于所确定的粒子数量变化情况调整并行分区边界位置。具体的，作为本发明的优选实例，所述步骤316包括：步骤3161、划定并行分区，并按照主方向自小至大的策略对每一并行分区所具有的粒子进行排序，以确定当前的并行分区所对应的粒子序列；步骤3162、在每一次时间步迭代循环后确定当前并行分区内的粒子数量变化情况：即自第一并行分区开始，判断当前区域内所对应的粒子是否多于初始粒子数，是则将粒子序列中多出的粒子划给下一并行分区(即第二并行分区)并更新下一并行分区的上边界和当前并行分区的下边界；步骤3163、重复步骤3162，以依次完成每一并行分区所对应的粒子数量变化情况，并基于所确定的粒子数量变化情况调整并行分区上下边界位置，如收到第一并行分区划分出的粒子后统计第二并行分区的粒子数，判断当前区域内所对应的粒子是否多于初始粒子数，是则将所确定的粒子序列中多出的粒子划给下一并行分区(即第三并行分区)并更新第三并行分区的上边界和第二并行分区的下边界，直至最后一个并行分区，此时仅统计所对应的粒子数量。

具体的，作为本发明的优选实例，所述步骤3还包括步骤32对每一并行分区进行区域搜索，其具体包括步骤321、将每一并行分区所对应的粒子，划分为该并行分区内部粒子以及分区边界通信粒子；划分规则为对当前行分区所对应的粒子进行遍历，并比较当前的粒子到该粒子所属并行分区的上边界或下边界的距离是否小于所设定的自身紧支域，若其到上边界的距离小于自身紧支域，则该粒子被确认为与上一并行分区进行通信的分区边界通信粒子；若其到下边界的距离小于自身紧支域，则该粒子被确认为与下一并行分区进行通信的分区边界通信粒子，若当前的粒子到该粒子所属并行分区的上边界、下边界的距离均不小于所设定的自身紧支域，则被确认为该并行分区内部粒子；步骤322、分别对每一并行分区的内部粒子进行搜索，对应的搜索策略为：步骤3221、确定当前并行分区所对应的计算域大小；步骤3222、按照一定比例将所确定的计算域划分为若干区域；步骤3223、将每一区域均划分为若干相同的背景网格单元，确定搜索方向并基于所确定的搜索方向计算各区域所对应的网格单元数，如果所述网格单元数大于所设定的上限阈值，则进一步缩小前述比例并重复步骤3222；步骤3224、确定搜索域并搜索当前区域所对应的临近粒子；所述搜索域的大小大于当前待搜索的区域大小，以使得所述搜索域能够与当前待搜索的区域的临近区域部分重叠；步骤3225、循环重复步骤3224直至搜索完成全部区域的搜索，并将搜索过程中所获得的各粒子信息进行更新；步骤3226、将当前搜索完成后的并行分区的下边界的边界粒子信息传递给下一并行分区。

需要说明的是：本部分在对每一并行分区或者称为计算节点的内部粒子进行搜索时，在这个节点上的搜索也是分区域(即步骤3222所述的区域)进行的，这个分区域搜索与前述并行分区的搜索不同，此时我们假设前计算节点内的粒子就是全部粒子，其通过设定搜索域或称为窗口搜索域搜索当前区域所对应的临近粒子，且每次只需要搜索窗口搜索域中的粒子，结束后依次移动窗口搜索域再搜索其他窗口中的粒子，这样做的好处就是不需要一次把计算节点内的所有粒子的包围盒都划上背景网格单元，每次只需要在窗口中有粒子的情况下进行搜索，而若窗口搜索域要比背景网格单元大很多，则在粒子分布极其不均匀的情况下大量没有粒子的窗口搜索域会一次跳过，节约大量的搜索时间，并且每一次的背景网格单元仅在窗口搜索域中划分，而窗口搜索域的大小又比整体包围盒小很多，则可以节约很多存储空间，这样在保证具有树搜索的部分优点的同时，又不需要面临递归过程中出现的栈溢出或递归过深的问题；实验证实，最佳的搜索域数量应该与每个搜索域中的背景网格的数量相当。同时在搜索域进行搜索，为了使得各搜索域之间能够搜索，需要使得搜索域各搜索域的每次搜索时都重复覆盖部分其他搜索域的内容，即步骤3224所述的确定搜索域并搜索当前区域所对应的临近粒子时，所述搜索域的大小大于当前待搜索的区域大小，以使得所述搜索域能够与当前待搜索的区域的临近区域部分重叠，如图3所示的例图，计算域被划分为均等的4块区域，在对任意一块区域进行搜索时，将每一区域均划分为若干层相同的背景网格单元且使得该搜索域略大于该区域即其与其他区域均存在一部分重叠交叉的区域，优选的可使得所述重叠交叉的区域为一层背景网格单元，如图4所示，其标注有横线的区域是当前待搜索的区域本身的大小，标注有斜线的区域是额外需要搜索的区域，以保证各区域间实现搜索；为了搜索域搜索不存在搜索重复的问题，我们借用有向搜索的方法，例如如图5所示例子，使得搜索沿着图5a所示的方向进行搜索时，如果当前搜索的是同一的背景网格单元(以中心的那个背景网格为例)，则需要判断该背景网格对应的节点号的大小，按照统一的大小标准判断(例如都是被搜索的节点号比较大的时候是一对作用对)避免重复搜索；当中心节点(即目标节点)位于左侧边界(搜索区域左侧大于当前分区的部分)的时候搜索路径如图5b所示，并且不搜索本身所在的网格，以防止重复搜索；当中心节点位于下方边界的时候，搜索路径如图5c，同样不搜索自身所在网格；当中心节点位于角落边界中，搜索路径如图5d，只有一个搜索网格，由上例可知同网格间的搜索不需要对比节点号就可以保证不重复搜索，这样做可以节省很多时间和计算资源，相比传统的链表搜索，速度不会减少，其内存的需求大量降低。

(4)、计算仿真模型所对应的状态方程以及控制方程，具体的，作为本发明的优选实例，在进行计算时，若仿真模型存在对称边界即所述仿真模型为对称模型，则对仿真模型进行作用对信息增补处理；所述进行作用对信息增补处理并不需要进行额外搜索，只需确认相应的镜像关系即可，具体的：

所述对仿真模型进行作用对信息增补处理包括：

步骤41、分别独立确认每一并行分区是否存在一个对称边界，若当前的并行分区存在对称边界，则对当前的并行分区中所对应的所有作用对进行循环，并判断某一作用对中是否存在至少一个粒子到对称边界的距离小于自身紧支域，是则确认可能存在具有镜像作用的粒子对；

步骤42、进一步确定该作用对的镜像关系，因为是对称的关系，所以只需要考虑一个粒子的镜像是否是另一个粒子的作用对，即判断该作用对中的任意一种粒子的镜像粒子是否是另一粒子的作用对，所述的判断策略为：获取该作用对中的任意一种粒子的镜像粒子的镜像位置，并判断所获取的镜像位置与另一粒子的间距与两个粒子紧支域之和的差值的绝对值是否小于所确定的判断阈值，是则确定该作用对为具有镜像作用的粒子对即确认该作用对中的任意一个粒子的镜像粒子与另一粒子互为作用对，且若该作用对中存在其到分区边界的距离小于自身紧支域，则确认该粒子的镜像粒子与该粒子互为作用对；依次对所有作用对循环，以获得全部的关于对称边界的作用对信息，对仿真模型进行作用对信息的进行增补。

进一步的，作为本发明的优选，由于所述对仿真模型不仅存在1/2对称，而且存在多重对称边界，如当对称边界为两个也就是说是1/4对称的时候，第一次产生的镜像位置对于第二个对称面同样会产生一个镜像位置，这个镜像位置同样有可能是一个作用粒子的位置，此外还有1/8对称的情况，针对上述三种情况的对称边界，可按照下述过程进行作用对信息增补处理；

步骤41′、分别独立确认每一并行分区是否存在对称边界，若当前的并行分区存在对称边界且对称边界个数至少为2个，则分别沿着x方向、y方向以及z方向对当前的并行分区中所对应的所有作用对进行循环，并在依次对上述3个方向进行循环的过程中，判断当前的作用对中是否存在至少一个粒子到其中任意一个对称边界的距离均小于自身紧支域，是则确认可能存在具有镜像作用的粒子对；

步骤42′、进一步确定该作用对的镜像关系，因为是对称的关系，所以需要考虑一个粒子在各个对称面上的镜像是否均是另一个粒子的作用对，即判断该作用对中的任意一种粒子的在各个对称面上的镜像粒子是否均是另一粒子的作用对，所述的判断策略为：首先获取该作用对中的任意一种粒子在某一对称面上镜像粒子的镜像位置，并判断所获取的镜像位置与另一粒子的间距与两个粒子紧支域之和的差值的绝对值是否小于所确定的判断阈值，是则确定确认该粒子在当前的对称面上镜像粒子与另一粒子互为作用对，且若该作用对中存在其到分区边界的距离小于自身紧支域，则确认该粒子在当前的对称面上的镜像粒子与该粒子互为作用对；其次获取该作用对中的任意一种粒子在其他对称面上镜像粒子的镜像位置，并判断所获取的镜像位置与另一粒子的间距与两个粒子紧支域之和的差值的绝对值是否小于所确定的判断阈值，是则确定确认该粒子在所述对称面上镜像粒子与另一粒子互为作用对，且若该作用对中存在其到分区边界的距离小于自身紧支域，则确认该粒子在所述对称面上的镜像粒子与该粒子也互为作用对；依次对所有作用对循环，以获得全部的关于对称边界的作用对信息，对仿真模型进行作用对信息的进行增补；同时所述作用对信息包括相互作用的粒子的标号以及各自所对应的作用方向，并在获取作用对信息之后，根据对称记号计算人工粘性和相互作用函数，同样记录在作用对信息之中。

上述对多重对称边界的处理过程可由下述步骤及公式表示：

(1)分别对x方向、y方向、z方向进行循环，确认是否存在对称边界；即

对z方向循环：for i＝0到s_z

对y方向循环：for j＝0到s_y

对x方向循环：for k＝0到s_x

(2)判断作用对中的一系列A点在各个对称面上的镜像粒子一系列A点是否均是点B(X_b，Y_b，Z_b)的作用对，其中一系列A点为坐标集合，可表示为(sign(k)*Xa，sign(j)*Ya，sign(i)*Za)；其中，x方向、y方向、z方向的对称记号分别为s_x，s_y，s_z，0表示没有对称边界，1表示有对称边界；sign为对称方向符号且sign(0)＝1，sign(1)＝-1；*表示相乘；是则记录A点以及相应的对称方向符号sign为一对作用对并计算这一对对应的核函数信息。

(3)结束x方向循环；结束y方向循环；结束z方向循环。

具体的，作为本发明的优选实例，所述步骤4还包括计算仿真模型所对应的物理控制方程，所述物理控制方程的计算包括针对固体相关部分计算以及流体相关部分计算，对应的，所述固体相关部分计算包括：a1、连续性方程计算，如首先对所有作用对循环，以获得作用粒子的坐标，速度和对称符号，通过对称符合合成计算用粒子的坐标和速度；然后获得这对作用对的相互作用函数dwdx；用合成后的相对速度乘以作用函数获得这对作用对的两个粒子的密度变化率增量；累加循环计算中每个粒子获得的增量后；结束循环，使得各并行分区间进行通信；b1、动量守恒方程计算：对所有作用对循环，获取各作用粒子的应力，对称符号，通过对称符号合成计算用粒子应力；然后获得各作用对相互作用的函数dwdx；用合成后的应力计算每一对粒子中第一个粒子的加速度增量，以这个增量为基础通过对称符号合成第二个粒子的加速度增量；使得每个粒子独自累加循环计算获得各自所对应的增量；结束循环，结束循环，使得各并行分区间进行通信并合成全部粒子的加速度；c1、能量守恒方程计算：对所有作用对循环，获取作用粒子的应力，速度，对称符号，通过对称符号合成计算用粒子应力和速度；获得各作用对相互作用的函数dwdx；用合成后的应力和速度计算每一对粒子中粒子的比内能变化率增量；使得每个粒子独自累加循环计算获得的增量；结束循环，使得各并行分区间进行通信并合成全部粒子的比内能变化率；对应的，所述流体相关部分计算包括：a2、动量守恒方程计算：对所有作用对循环，获取作用粒子的压力和这个作用对的人工粘性，对称符号，当不采用人工粘性时，通过两个粒子的压力合成黎曼解的压力；获得各作用对相互作用的函数dwdx；依次计算每一对粒子中第一个粒子的加速度增量，以这个增量为基础通过对称符号合成第二个粒子的加速度增量；使得每个粒子独自累加循环计算获得的增量；结束循环，使得各并行分区间进行通信并合成全部粒子的加速度；b2、能量守恒方程计算：对所有作用对循环，获取每一作用粒子的压力和这个作用对的人工粘性，速度，对称符号，通过对称符号合成计算用粒子速度，当不采用人工粘性时，通过压力合成黎曼解的压力；获得各作用对相互作用的函数dwdx；用合成后的速度计算每一对粒子的比内能变化率增量；使得每个粒子独自累加循环计算获得的增量；结束循环，使得各并行分区间进行通信并合成全部粒子的比内能变化率。

具体的，作为本发明的优选实例，所述步骤4还包括通过进行时间积分以更新速度和坐标以及其他状态和物理量：如对速度和坐标的更新：对所有并行分区的所有实粒子进行循环；获取各粒子的加速度、当前速度和坐标；使得各粒子更新后的速度等于当前速度加上加速度与时间步长的乘积，更新后的坐标等于当前坐标加上更新后的速度与时间步长的乘积；如对能量变化的更新：对所有并行分区的所有实粒子循环，获取各粒子的当前比内能和比内能变化率，使得更新后各粒子的比内能等于当前比内能加上比内能变化率与时间步长的乘积；如对应变的更新：对所有并行分区的所有实粒子进行循环；获取各粒子当前的应变以及应变率；对每个应变分量循环，使得更新后的应变分量等于当前应变分量加上对应的应变率与时间步长的乘积；如对失效状态的更新：对所有并行分区的所有实粒子进行循环，获取粒子更新后的应力和应变代入所设定的公式判断粒子失效状态，如果失效则标记为失效粒子；如对侵蚀状态的更新：对所有并行分区的所有实粒子进行循环，获取各粒子的时间步长，应变等信息，代入对应的公式判断侵蚀状态，如果达到侵蚀标准则标记为侵蚀粒子，再次对所有粒子循环并删除掉有侵蚀标记的粒子，最后将未侵蚀的粒子重新放入连续的存储数据结构中作为下一步计算的粒子数据。

具体的，作为本发明的优选实例，所述步骤4还包括进行边界条件计算即根据所设定的边界条件，处理受到约束的粒子的速度和进行坐标更新，具体包括下述步骤：对所有并行分区的所有实粒子循环；判断当前粒子是否是被标记为边界标记的粒子，如果粒子不是边界标记的粒子，则正常更新该粒子的速度和坐标；如果粒子是边界标记的粒子，则进一步判断该粒子的位移是否达到边界条件的约束，如果达到约束，则该粒子的位移和坐标不再更新，且速度置为0如果没有达到约束，则该粒子的速度和坐标正常更新。

具体的，以爆轰相关的传爆路径搜索为例进行说明：若采用直接搜索方法：

则首先确定该模型的初始设置的起爆位置；然后获取起爆时间和当前时间的时间差；获取炸药的爆速并确定起爆条件，所述起爆条件的确定是指如果粒子到起爆点位置小于或等于爆速乘以时间差，则达到起爆条件，否则没有达到起爆条件；最后遍历所有粒子，如果当前粒子达到起爆条件，则设置该粒子为起爆粒子，否则设置为未起爆粒子。若采用间接搜索方法：由于间接搜索不需要初始起爆位置，其是通过粒子与已起爆粒子的相对位置关系确定，因此其包括：对该模型的除了镜像作用对之外的所有作用对循环，依次判断是否一对作用对中有一个已经起爆粒子和一个未起爆粒子，如果不满足前述条件则继续判断下一个作用对；如果满足条件则获取其中一个起爆粒子的起爆时间，计算与当前时间的时间差，用时间差乘以爆速获取爆轰波传播的距离，并进一步判断这个距离与这对作用对之间的距离，如果粒子间距离大于爆轰波传播距离，则进行下一对作用对判断，否则未起爆粒子将被起爆，记录临时起爆记号(临时起爆记号不会影响后面搜索的判断)；最后完成所有作用对循环之后采用临时起爆记号更新粒子起爆记号，使得各并行分区进行通信，至此完成起爆传爆的搜索过程。

(5)、对所述仿真模型进行后处理，即基于所设定终止条件，输出对应的仿真结果。具体的，作为本发明的优选实例，所述步骤5包括：若当前计算的时间满足输出条件的时候，将所有并行分区所对应的仿真数据合并到所设定的主进程下并获取当前输出文件的序号标记；对前述主进程中所有粒子进行循环；将各粒子对应输出的变量及变量名称写入文件；完成显示需要的数据之后，将重启动需要的增量信息以及材料和计算控制信息写入同一文件；最后输出对应的仿真结果，且更新当前文件的输出序号标记。具体的，作为本发明的优选实例，所述步骤5包括设置屏幕打印计算进度，即：对所有并行分区进行通信，获得各个并行分区所对应的屏幕打印需要参数的值；对各所获得的值进行筛选，获取其中的最大值，最小值以及各值的平均值；将筛选后的各值代入所设定的屏幕打印的函数；根据所设定的输出的行数调整光标位置，令每一次的刷新刚好覆盖需要刷新的内容，以防止出现错行的问题；执行屏幕打印任务。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。

Claims (1)

1.一种用于SPH算法的对称边界处理方法，其特征在于，包括如下步骤：

步骤1、创建待进行仿真分析的仿真模型并对所创建的仿真模型进行初始化设置；

所述的步骤1其包括：

步骤101、创建仿真模型并设定所创建的仿真模型所对应的网格文件，获得对应的网格信息并存入临时的网格数据中，所述网格信息包括网格序号以及网格材料号，且网格信息包含与其所对应的节点号；

步骤102、设定所创建的仿真模型所对应的节点文件，获得对应的节点信息并存入临时的节点数据中，所述节点信息包括节点号以及节点坐标；

步骤103、结合网格信息与节点信息计算仿真模型所对应的各个网格的中心坐标以及每一网格所对应的体积数据，并基于所计算的体积数据计算各自对应的等效直径；

步骤104、将各网格转换为所对应的粒子，其中每一粒子的位置即其所对应的网格的中心坐标，其大小即其所对应的网格的等效直径，其材料号即其所对应的网格材料号；

步骤105、设定材料文件，获取仿真模型所对应的每种材料的密度，泊松比，弹性模量以及材料模型的类型和参数；

步骤106、根据每一粒子的所对应的网格材料号将所设定的材料的属性和模型参数赋给每一粒子；

步骤107、设定边界条件和初始条件，获取仿真模型所对应的条件值和条件作用的组件号，存在临时的初边值数据中；

步骤108、设定组件信息文件并获得仿真模型所对应的组件数，若按顺序打开组件文件，则依次设定对应的组件号、组件类型；进一步的，若当前的组件号和组件类型满足当前所对应的初边界条件，则继续设定当前组件的成员，获取相关网格的序号，将初边值条件赋给网格相应的粒子；若边界条件作用于节点，则获取该节点属于的单元数据，将边界条件赋给对应单元转换得到的粒子；

步骤109、设定初始控制参数；所述初始控制参数包括当前时间、总体计算时长、当前时间步长、总体循环数、当前时间步数、最小时间步长、最大时间步长、对称信息和保存结果设定、并行信息；

步骤1010、设定每种材料的属性和模型参数，将所设定的材料的属性和模型参数按照粒子的材料号赋给每个粒子；

步骤2、计算经初始化设置后的仿真模型所对应的每一粒子的压力值、应力张量以及人工粘性值；

步骤3、在所设定的每一个时间步内对所述仿真模型进行搜索分析，以确定每一粒子所对应的近邻粒子；所述的搜索分析包括：

步骤31确定仿真模型所对应的各并行分区以及各并行分区边界，其包括：

步骤311、首先对仿真模型所对应的所有粒子进行循环，分别找到在三维空间的X方向、Y方向、Z方向上，坐标最小所对应的粒子I以及坐标最大所对应的粒子II；其次将粒子II的坐标与粒子I的坐标相减，得到X、Y、Z三个方向上的最大跨度值；最后通过将所述仿真模型所对应的每一粒子的坐标依次与粒子I的坐标相减并除以所获得的该坐标方向上的最大跨度值，确定出每一粒子的相对无量纲坐标(X_j-X_i/X₀,Y_j-Y_i/Y₀,Z_j-Z_i/Z₀),其中，X_j、Y_j、Z_j分别表示X方向、Y方向、Z方向上的粒子坐标；X_i、Y_i、Z_i分别表示X方向、Y方向、Z方向上坐标最小的粒子I的坐标，X₀表示X方向上的最大跨度值，Y₀表示Y方向上的最大跨度值，Z₀表示Z方向上的最大跨度值，j≥1；

步骤312、分别计算所有粒子所对应的无量纲坐标在X方向、Y方向、Z方向上的平方和，并将平方和最小的无量纲坐标所构成的向量方向确定为主方向；

步骤313、基于所确定的主方向，按照所有粒子与主方向的相对坐标从小至大进行排序；

步骤314、确定待划定的并行分区数量并确定每一并行分区所对应的初始粒子数以及各并行分区边界位置所对应的分界粒子序号；

步骤315、确定每一并行分区边界位置；

步骤316、划定并行分区并对每一并行分区的粒子进行排序，并通过各并行分区通信来确定各个并行分区内的粒子数量变化情况，并基于所确定的粒子数量变化情况调整并行分区边界位置；

步骤4、计算仿真模型所对应的状态方程以及控制方程，在进行计算时，若仿真模型存在对称边界，则对仿真模型进行作用对信息增补处理；

若所述对仿真模型为1/2对称模型，所述对仿真模型进行作用对信息增补处理包括：

步骤41、分别确认每一分区是否存在对称边界，若当前的分区存在对称边界，则对当前的分区中所对应的所有作用对进行循环，并判断某一作用对中是否存在至少一个粒子到对称边界的距离小于自身紧支域，是则确认可能存在具有镜像作用的粒子对；

步骤42、进一步确定该作用对的镜像关系，即判断该作用对中的任意一种粒子的镜像粒子是否是另一粒子的作用对，判断策略为：获取该作用对中的任意一种粒子的镜像粒子的镜像位置，并判断所获取的镜像位置与另一粒子的间距与两个粒子紧支域之和的差值的绝对值是否小于所确定的判断阈值，是则确定该作用对为具有镜像作用的粒子对即确认该作用对中的任意一个粒子的镜像粒子与另一粒子互为作用对，且若步骤41所述的某一作用对中存在其到分区边界的距离小于自身紧支域，则确认该粒子的镜像粒子与该粒子互为作用对；依次对所有作用对循环，以获得全部的关于对称边界的作用对信息，对仿真模型进行作用对信息增补处理；

若所述对仿真模型为多重对称模型即1/4对称模型或者1/8对称模型，所述对仿真模型进行作用对信息增补处理包括：

步骤41′、分别独立确认每一并行分区是否存在对称边界，若当前的并行分区存在对称边界且对称边界个数至少为2个，则分别沿着x方向、y方向以及z方向对当前的并行分区中所对应的所有作用对进行循环，并在依次对x方向、y方向以及z方向进行循环的过程中，判断当前的作用对中是否存在至少一个粒子到其中任意一个对称边界的距离均小于自身紧支域，是则确认可能存在具有镜像作用的粒子对；

步骤42′、进一步确定该作用对的镜像关系，因为是对称的关系，所以需要考虑一个粒子在各个对称面上的镜像是否均是另一个粒子的作用对，即判断该作用对中的任意一种粒子的在各个对称面上的镜像粒子是否均是另一粒子的作用对，判断策略为：首先获取该作用对中的任意一种粒子在某一对称面上镜像粒子的镜像位置，并判断所获取的镜像位置与另一粒子的间距与两个粒子紧支域之和的差值的绝对值是否小于所确定的判断阈值，是则确定确认该粒子在当前的对称面上镜像粒子与另一粒子互为作用对，且若该作用对中存在其到分区边界的距离小于自身紧支域，则确认该粒子在当前的对称面上的镜像粒子与该粒子互为作用对；其次获取该作用对中的任意一种粒子在其他对称面上镜像粒子的镜像位置，并判断所获取的镜像位置与另一粒子的间距与两个粒子紧支域之和的差值的绝对值是否小于所确定的判断阈值，是则确定确认该粒子在所述对称面上镜像粒子与另一粒子互为作用对，且若步骤41′所述当前的作用对中存在其到分区边界的距离小于自身紧支域，则确认该粒子在所述对称面上的镜像粒子与该粒子也互为作用对；依次对所有作用对循环，以获得全部的关于对称边界的作用对信息，对仿真模型进行作用对信息的进行增补；同时所述作用对信息包括相互作用的粒子的标号以及各自所对应的作用方向，并在获取作用对信息之后，根据对称记号计算人工粘性和相互作用函数，同样记录在作用对信息之中；

步骤5、基于所设定终止条件，输出对应的仿真结果。