一种基于粒子分裂的sph计算优化方法及系统
技术领域
1.本发明属于粒子优化领域，尤其涉及一种基于粒子分裂的sph计算优化方法及系统。
背景技术：
2.近年来，随着sph方法的飞速发展，引入“粒子分裂和聚合”成功实现计算区域内sph粒子的“空间分辨率”可变，使得计算一些大尺度规模的流动问题成为了可能。变分辨率技术是提高sph计算效率的一种有效方法，它可以在重点关注的核心计算区域使用高分辨率，而在其它区域使用低分辨率。变分辨率技术可以有效的平衡计算效率及计算效果，在传统网格方法中已获得了较为广泛的应用。对于sph这种纯lagrange粒子法，粒子优化是改变计算精度的有效途径。现有的sph方法中没有出现在进行粒子优化时，采用粒子分裂进行高效处理和计算的方法。关于粒子优化处理主要有粒子重置、自适应增减粒子等。
3.但是，总体上来看，两类粒子优化算法在计算精度、稳定性及计算效率等方面还存在如下诸多问题。第一，两类算法均缺乏有效的手段以保证大-小粒子作用区域的计算精度和稳定性，现有粒子优化算法中大-小粒子通过分段连续的复合核函数直接相互作用，这些核函数在不同区间使用不同的光滑长度，因此，核函数本身并不光滑，sph算法的插值精度和计算稳定性非常低。第二，现有粒子优化算法的粒子优化参数选取较为困难，粒子优化参数包括分离系数ε及光滑比α，前者决定着大-小粒子之间的分离程度，后者决定着大-小粒子的光滑长度之比。现有算法通过理论分析或算例测试给出ε与α的取值建议普适性较差，对于某些数值问题，同时确定“理想的”、能达到最优效果的ε与α的值非常困难。
技术实现要素：
4.本发明旨在解决上述问题，提供一种基于粒子分裂的sph计算优化方法及系统。
5.第一方面，本发明提供一种基于粒子分裂的sph计算优化方法，包括大粒子、小粒子和代理粒子；包括如下步骤：所述代理粒子指在优化过程中发生作用的大粒子与小粒子；所述代理粒子与在优化过程中和小粒子发生所作用的大粒子相关联；定义代理粒子是为保证粒子优化算法的计算精度，防止由于大-小粒子的直接作用产生计算误差；步骤一、将符合粒子优化准则的大粒子分裂为小粒子；所述优化准则一种为划定优化区域，一种为选取粒子类型；两种优化准则可根据实际优化需求选择一种或同时采用两种。
6.步骤二、重置代理粒子的物理参数；包括所述代理粒子的质量、粒子间距与前述大粒子分裂时产生的小粒子相同；在数值模拟过程中，当大粒子与小粒子发生作用时，大粒子的关联代理粒子将取代大粒子并与小粒子相互作用，将所述代理粒子的密度、速度、压力、光滑长度与代理粒子相关联的大粒子相同；步骤三、遍历计算域搜索粒子对，并确定粒子在当前时间步所属的计算域；
步骤四、通过求解连续性方程、动量方程及能量方程，获得大粒子、代理粒子及小粒子所承载的粒子密度dρ、粒子速度dv；步骤五、通过下式将代理粒子所承载的粒子密度dρ、粒子速度dv参量进行平均并传递给与代理粒子相关联的大粒子：式中，、分别代表过渡区大粒子及其关联的代理粒子的物理参量；n为粒子分裂的数量；一般本发明所述计算优化方法中n=ke，其中k≥2；e取值2或3，在二维情况下取值为2，三维情况下取值为3；步骤六、更新大粒子、小粒子的粒子参数，当前时间步计算结束；判断时间步次数是否大于预设次数；若否，则继续重复前述步骤进行下一时间步计算；若是，则优化结束。
7.进一步，本发明所述基于粒子分裂的sph计算优化方法，所述计算域包括三个子区域，分别为大粒子区、过渡区及小粒子区；所述大粒子区为不与小粒子发生作用的大粒子构成的区域；所述过渡区为与小粒子发生作用的大粒子构成的区域；所述小粒子区为由小粒子构成的区域。
8.进一步，本发明所述基于粒子分裂的sph计算优化方法，其特征在于：所述遍历计算域搜索粒子对时，按不同子区域搜索相应粒子对，并确定当前时间步粒子所属的子区域。粒子属于子区域的信息将用于下一时间步中的粒子对搜索及粒子信息计算。
9.进一步，本发明所述基于粒子分裂的sph计算优化方法，为保证在粒子分裂过程中的质量、动量及动能守恒，所述大粒子和小粒子的质量m、密度ρ及速度v设定为：密度ρ及速度v设定为：密度ρ及速度v设定为：；式中，mf为小粒子质量、mc为大粒子质量、ρf为小粒子密度、ρc为大粒子密度、vf为小粒子速度、vc为大粒子速度；所述大粒子与分裂形成的小粒子之间的间距、及小粒子的光滑长度hf设定为：定为：式中，dpc表示大粒子的初始粒子间距；hc为大粒子的光滑长度；α为光滑比，定义大-小粒子的光滑长度之比；ε为分离系数，定义大-小粒子间的水平和竖直距离。
10.进一步，本发明所述基于粒子分裂的sph计算优化方法，所述大粒子区内粒子的临近粒子是大粒子；大粒子区域内某一粒子i的临近粒子j位于过渡区时，j粒子以大粒子形式
参与i粒子相关物理参数的计算；所述过渡区内代理粒子的临近粒子是小粒子或大粒子代理粒子；过渡区外大粒子j位于过渡区内代理粒子i的支持域内时，j粒子的代理粒子参与计算代理粒子i的相关物理参数；所述小粒子区的小粒子为建模时指定的小粒子或由大粒子分裂产生的小粒子；小粒子的临近粒子是小粒子或过渡区的代理粒子。
11.进一步，本发明所述基于粒子分裂的sph计算优化方法，所述大粒子区内的流场信息由光滑长度为hc的核函数插值计算；所述过渡区内的流场信息通过大粒子的代理粒子计算，由光滑长度为hr=hc的核函数插值计算；所述小粒子区内的流场信息由光滑长度为hf=αhc的核函数插值计算。
12.进一步，本发明所述基于粒子分裂的sph计算优化方法，所述粒子分裂采用沿坐标轴向的正方形粒子分裂模型；所述正方形粒子分裂模型将大粒子分裂为四个小粒子，小粒子分布在以大粒子为中心的正方形的四个角上，正方形的边与坐标轴平行。
13.进一步，本发明所述基于粒子分裂的sph计算优化方法，所述正方形粒子分裂模型的分离系数ε为0.5。首先，ε=0.5有助于减小粒子分裂后小粒子之间的重叠，改善粒子秩序；其次，ε取固定值后，在数值模拟过程中，唯一需要调整的参数只有光滑比α，这大大提高了本发明所述粒子优化算法的实用性。
14.第二方面，本发明提供一种基于粒子分裂的sph计算优化系统，包分裂单元、参数重置单元、遍历搜索单元和分析输出单元；所述分裂单元用于将符合粒子优化准则的大粒子分裂为小粒子；所述参数重置单元用于重置代理粒子的物理参数；包括所述代理粒子的质量、粒子间距与前述大粒子分裂时产生的小粒子相同；所述代理粒子的密度、速度、压力、光滑长度与代理粒子相关联的大粒子相同；所述遍历搜索单元用于遍历计算域搜索粒子对，并确定粒子在当前时间步所属的计算域；所述分析输出单元用于通过求解连续性方程、动量方程及能量方程，获得大粒子、代理粒子及小粒子所承载的粒子密度dρ、粒子速度dv；并通过下式将代理粒子所承载的粒子密度dρ、粒子速度dv参量进行平均并传递给与代理粒子相关联的大粒子：；式中，、分别代表过渡区大粒子及其关联的代理粒子的物理参量；然后更新大粒子、小粒子的粒子参数，当前时间步计算结束；判断时间步次数是否大于预设次数；若否，则继续重复前述步骤进行下一时间步计算；若是，则优化结束；输出优化结果。
15.第三方面，本发明提供一种基于粒子分裂的sph计算优化系统，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其特征在于：所述处理器执行所述程序时实现第一方面所述优化方法的步骤。
16.本发明所述基于粒子分裂的sph计算优化方法及系统，根据粒子优化准则，将满足
优化准则的大粒子分裂为小粒子，将大-小粒子的接触区域定义为过渡域，引入代理粒子的概念，利用代理粒子代替大粒子与小粒子进行作用。保持了sph方法在大-小粒子的接触区域的空间插值精度，避免了由于大-小粒子的直接接触而造成的插值误差和数值震荡，在保证计算精度的同时有效提高了数值模拟效率。
附图说明
17.图1为本发明实施例所述基于粒子分裂的sph计算优化方法粒子配置示意图；图2为本发明实施例所述基于粒子分裂的sph计算优化方法计算域示意图；图3为本发明实施例所述基于粒子分裂的sph计算优化方法流程示意图；图4为本发明实施例三所述基于粒子分裂的sph计算优化系统结构示意图；图5为本发明实施例四所述基于粒子分裂的sph计算优化系统结构示意图；图6为本发明实施例二所述小球落水算例粒子配置示意图；图7为本发明实施例二所述小球入水深度的仿真与实验结果对比示意图；图8为本发明实施例二所述不同粒子配置时流场形态及压力场随时间变化过程示意图。
具体实施方式
18.下面通过附图及实施例对本发明所述基于粒子分裂的sph计算优化方法及系统进行详细说明。
19.实施例一本公开实施例采用沿坐标轴向的正方形粒子分裂模型。该粒子分裂模型将大粒子分裂为四个小粒子，小粒子分布在以大粒子为中心的正方形的四个角上，正方形的边与坐标轴平行。为保证在粒子分裂过程中的质量、动量及动能守恒，对分裂后形成小粒子的质量、密度及速度定义如下：
ꢀꢀ
（1）
ꢀꢀ
（2）
ꢀꢀ
（3）式中，mf为小粒子质量、mc为大粒子质量、ρf为小粒子密度、ρc为大粒子密度、vf为小粒子速度、vc为大粒子速度；大粒子与分裂形成的小粒子之间的间距、及小粒子的光滑长度hf设定为：
ꢀꢀ
（4）
ꢀꢀ
（5）式中，dpc表示大粒子的初始粒子间距，hc为大粒子的光滑长度，α、ε分别称为光滑比和分离系数，分别定义大-小粒子的光滑长度之比和大-小粒子间的水平和竖直距离。
20.在本公开实施例中粒子分裂模型的分离系数ε固定为0.5，这样做有以下两方面的优势：首先，ε=0.5有助于减小粒子分裂后小粒子之间的重叠，改善粒子秩序；其次，ε取固定
值后，在数值模拟过程中，唯一需要调整的参数只有光滑比α，这大大提高本实施例所述粒子优化算法实用性及效率。
21.当粒子优化算法被引入sph计算框架后，sph方法的连续性会下降。因此，为保证粒子优化算法的计算精度，防止由于大-小粒子的直接作用产生计算误差，本发明所述基于粒子分裂的sph优化算法中引入了一种新型粒子“代理粒子”，并基于代理粒子提出了一种大-小粒子的混合作用格式。
22.在本公开实施例所述基于粒子分裂的sph计算优化方法中共包括三种类型的粒子：大粒子、小粒子及代理粒子。三种类型的粒子配置如图1所示，代理粒子是与大粒子相关联的，在本公开实施例中每个大粒子都关联着四个代理粒子。因此，在本公开实施例中的代理粒子的质量及粒子间距与粒子分裂时产生的小粒子相同，即：
ꢀꢀ
（6）
ꢀꢀ
（7）式中，下标r表示代理粒子的相关物理量，mr是代理粒子的质量，dpr、dpf分别为代理粒子与小粒子的初始粒子间距。
23.在数值模拟计算过程中，当大粒子与小粒子发生作用时，大粒子的关联代理粒子将取代大粒子并与小粒子相互作用，因此规定代理粒子的密度、速度、压力、光滑长度等参数与其关联的大粒子相同，即：
ꢀꢀ
（8）
ꢀꢀ
（9）
ꢀꢀ
（10）
ꢀꢀ
（11）本公开实施例所述sph计算优化方法中，整个计算域可以划分为三个子区域，分别为：大粒子区、过渡区及小粒子区。如图2所示，（a）为大粒子区，（b）为过渡区，（c）为小粒子区，三个子区域的定义及相互作用方式如下：（a）大粒子区：由不与小粒子发生作用的大粒子构成的区域。如图2所示，该区域内粒子的临近粒子只能是大粒子。当该区域内某一粒子i的临近粒子j位于过渡区时，j粒子以大粒子形式参与i粒子相关物理量的计算。该区域内的流场信息由光滑长度为hc的核函数插值计算。
24.（b）过渡区：由与小粒子发生作用的大粒子构成的区域。如图2所示，过渡区内的流场信息通过大粒子的代理粒子计算。代理粒子的临近粒子只能是小粒子或大粒子代理粒子。当过渡区外大粒子j位于过渡区内代理粒子i的支持域内时，j粒子的代理粒子参与计算代理粒子i的相关物理量。该区域内的流场信息由光滑长度为hr=hc的核函数插值计算。
25.（c）小粒子区：由小粒子构成的区域，这些小粒子可以是建模时指定的，也可以是模拟过程中由大粒子分裂产生的。如图2所示，小粒子的临近粒子只能是小粒子或过渡区的代理粒子。该区域内的流场信息由光滑长度为hf=αhc的核函数插值计算。
26.由以上分析可知，本公开实施例所述基于粒子分裂的sph计算优化方法中，过渡区的粒子可被看作混合粒子：当与大粒子区粒子相互作用时，该区域内的粒子是大粒子；当与小粒子发生作用或过渡区内粒子自身相互作用时，过渡区内的代理粒子将代替大粒子参与流场计算；这种处理方法防止了由于大-小粒子直接相互作用而带来的计算误差和数值不稳定。因此，我们将以上不同子区域内粒子相互作用方式也称为混合粒子作用格式。
27.本公开实施例所述基于粒子分裂的sph计算优化方法如图3所示，包括如下步骤：步骤一、将符合粒子优化准则的大粒子分裂为小粒子；本实施例中采取的优化准则为划定优化区域；步骤二、重置代理粒子的物理参数；包括所述代理粒子的质量、粒子间距与前述大粒子分裂时产生的小粒子相同；在数值模拟过程中，当大粒子与小粒子发生作用时，大粒子的关联代理粒子将取代大粒子并与小粒子相互作用，将所述代理粒子的密度、速度、压力、光滑长度与代理粒子相关联的大粒子相同；步骤三、遍历计算域搜索粒子对，并确定粒子在当前时间步所属的计算域；遍历计算域搜索粒子对时，按不同子区域搜索相应粒子对，并确定当前时间步粒子所属的子区域；粒子属于子区域的信息将用于下一时间步中的粒子对搜索及粒子信息计算；步骤四、通过求解连续性方程、动量方程及能量方程，获得大粒子、代理粒子及小粒子所承载的粒子密度dρ、粒子速度dv；通过下式将代理粒子所承载的粒子密度dρ、粒子速度dv参量进行平均并传递给与代理粒子相关联的大粒子：
ꢀꢀ
（12）式中，、分别代表过渡区大粒子及其关联的代理粒子的物理参量；步骤五、更新大粒子、小粒子的粒子参数，当前时间步计算结束；步骤六、判断时间步次数是否大于预设次数；若否，则继续重复前述步骤进行下一时间步计算；若是，则优化结束。
28.本公开实施例所述sph计算优化方法从sph方法的lagrange粒子本质出发，通过对某一关键区域进行粒子优化的手段实现变分辨率数值模拟，保持了sph方法在大-小粒子的接触区域的空间插值精度，避免了由于大-小粒子的直接接触而造成的插值误差和数值震荡，在保证计算精度的同时有效提高了数值模拟效率。
29.实施例二本实施例基于前述实施例一所述基于粒子分裂的sph计算优化方法，用水球落水算例测试本发明所述计算优化方法对于包含流体与刚体相互作用及自由表面流动问题的计算效果。
30.小球落水算例的模型及粒子配置如图6所示，其中水的宽度和深度分别为2m、0.5m，小球直径r=5.5
×
10-2
m，重力加速度g=-9.81m/s2。计算开始时，小球刚接触水面，根据grennhow的实验数据，此时小球速度为v0=2.955m/s，声速cs=10v0=29.55m/s，人工粘性参数α
av
=0.05。小球落水过程的关键是小球与水之间接触力的计算，因此，本实施例中以小球为中心，设置边长为6r的粒子优化域；在计算过程中，优化域随着小球的运动而移动。水采用
40000个大粒子离散，粒子间距dpc=5.0
×
10-3
m，光滑长度hc=1.1dpc。
31.本公开实施例分别对比了以下粒子配置时的计算效果：本发明所述基于粒子分裂的sph优化算法，ε=0.5，α=0.5、0.6、0.7、0.75；reyes l
ó
pez粒子优化算法，ε=0.5，α=0.7；vacondio粒子优化算法，ε=0.4，α=0.9；传统单精度sph算法，流体粒子数n=4096、10000、40000、160000。
32.如图7所示为小球入水深度的计算结果与实验结果的对比，传统单精度sph算法的计算结果随着粒子数的增大而逐渐接近实验值。α=0.7、0.75时本发明所述粒子优化算法的计算结果与n=160000时单精度sph方法的计算结果非常接近。由于优化域随着小球的运动而移动，计算过程中的粒子数随着粒子优化而不断增加，因此，本实施例中对比了不同粒子配置情况下，计算50000个时间步的cpu总耗时（见表1）。
33.表 1 小球落水算例：不同粒子配置的计算效率对比由表1可见，本实施例与reyes l
ó
pez粒子优化算法的流体粒子数均在60000以下，与n=160000时单精度sph方法相比，本实施例与reyes l
ó
pez粒子优化算法节省了60%以上的计算时间。
34.为进一步验证粒子优化算法的计算效果，本实施例还对仿真获得的小球落水过程中的流场形态进行了分析。图8为传统单精度sph方法在流体粒子数n=160000时计算出的流场形态与压力与本发明所述方法在n=40000、α=0.7时的计算结果对比，（其中，每幅图像的左边为传统单精度sph方法使用160000个流体粒子时的仿真结果；右边为使用40000个流体粒子和新型粒子优化算法α=0.7时的仿真结果，虚线标出的是粒子优化域）可以看出，在小球落水的各个时刻，本发明所述方法与传统sph方法的计算结果在流场形态、压力波形态等方面非常一致，本发明所述方法计算出的压力场在大-小粒子区域之间过渡平滑。挡光滑比α变化时，压力场与速度场的变化一致，α=0.5的计算结果在大-小粒子接触区域存在轻微震荡，α=0.6、0.75时的计算结果与α=0.7时的计算结果基本一致。
35.由以上算例分析可见，对于小球落水问题，α=0.7时本发明所述方法的计算效率略低于reyes l
ó
pez方法，但是，α=0.6、0.7、0.75时，本发明所述方法在小球落水深度及压力场、速度场等方面的计算结果更为准确。
36.实施例三基于与前述实施例中一种基于粒子分裂的sph计算优化同样发明构思，本实施例
公开了一种基于粒子分裂的sph计算优化系统，如图4所示，所述优化系统包分裂单元、参数重置单元、遍历搜索单元和分析输出单元；其中所述分裂单元用于将符合粒子优化准则的大粒子分裂为小粒子；所述参数重置单元用于重置代理粒子的物理参数；包括所述代理粒子的质量、粒子间距与前述大粒子分裂时产生的小粒子相同；所述代理粒子的密度、速度、压力、光滑长度与代理粒子相关联的大粒子相同；所述遍历搜索单元用于遍历计算域搜索粒子对，并确定粒子在当前时间步所属的计算域；所述分析输出单元用于通过求解连续性方程、动量方程及能量方程，获得大粒子、代理粒子及小粒子所承载的粒子密度dρ、粒子速度dv；并通过下式（12）将代理粒子所承载的粒子密度dρ、粒子速度dv参量进行平均并传递给与代理粒子相关联的大粒子： （12）；式中，、分别代表过渡区大粒子及其关联的代理粒子的物理参量；然后更新大粒子、小粒子的粒子参数，当前时间步计算结束；判断时间步次数是否大于预设次数；若否，则继续重复前述步骤进行下一时间步计算；若是，则优化结束；输出优化结果。
37.通过对前述基于粒子分裂的sph计算优化方法的详细描述，本领域技术人员可以清楚的知道本实施例所述基于粒子分裂的sph计算优化系统的实施方法，为了说明书的简洁，在此不再详述。
38.实施例四基于与前述实施例中基于粒子分裂的sph计算优化方法的发明构思，本发明还提供一种基于粒子分裂的sph计算优化系统，包括处理器、存储器，所述存储器上存储有计算机程序，该程序被处理器执行时实现实施例一所述基于粒子分裂的sph计算优化方法的任意步骤。
39.图5为本实施例所述基于粒子分裂的sph计算优化系统的结构示意图。其中，总线可以包括任意数量的互联的总线和桥，总线将包括由处理器代表的一个或多个处理器和存储器代表的存储器的各种电路链接在一起。总线还可以将诸如外围设备、稳压器和功率管理电路等之类的各种其他电路链接在一起，这些都是本领域所公知的，因此，此处不再对其进行进一步描述。总线接口在总线和接收器和发送器之间提供接口。接收器和发送器可以是同一个元件，用于在传输介质上与各种其他系统通信的单元。处理器负责管理总线和通常的处理，而存储器被用于存储处理器在执行操作时所使用的数据。
40.本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质（包括但不限于磁盘存储器、cd-rom、光学存储器等）上实施的计算机程序产品的形式。
41.本发明是参照根据本发明实施例的方法、设备（系统）、和计算机程序产品的流程
图来描述的。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程中指定的功能的系统。
42.这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令系统的制造品，该指令系统实现在流程图一个流程或多个流程一个方框或多个方框中指定的功能。
43.这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程中指定的功能的步骤。
44.尽管已描述了本发明的优选实施例，但本领域内的技术人员一旦得知了基本创造性概念，则可对这些实施例进行另外的变更和修改，若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。
技术领域
1.本发明属于粒子优化领域，尤其涉及一种基于粒子分裂的sph计算优化方法及系统。
背景技术：
2.近年来，随着sph方法的飞速发展，引入“粒子分裂和聚合”成功实现计算区域内sph粒子的“空间分辨率”可变，使得计算一些大尺度规模的流动问题成为了可能。变分辨率技术是提高sph计算效率的一种有效方法，它可以在重点关注的核心计算区域使用高分辨率，而在其它区域使用低分辨率。变分辨率技术可以有效的平衡计算效率及计算效果，在传统网格方法中已获得了较为广泛的应用。对于sph这种纯lagrange粒子法，粒子优化是改变计算精度的有效途径。现有的sph方法中没有出现在进行粒子优化时，采用粒子分裂进行高效处理和计算的方法。关于粒子优化处理主要有粒子重置、自适应增减粒子等。
3.但是，总体上来看，两类粒子优化算法在计算精度、稳定性及计算效率等方面还存在如下诸多问题。第一，两类算法均缺乏有效的手段以保证大-小粒子作用区域的计算精度和稳定性，现有粒子优化算法中大-小粒子通过分段连续的复合核函数直接相互作用，这些核函数在不同区间使用不同的光滑长度，因此，核函数本身并不光滑，sph算法的插值精度和计算稳定性非常低。第二，现有粒子优化算法的粒子优化参数选取较为困难，粒子优化参数包括分离系数ε及光滑比α，前者决定着大-小粒子之间的分离程度，后者决定着大-小粒子的光滑长度之比。现有算法通过理论分析或算例测试给出ε与α的取值建议普适性较差，对于某些数值问题，同时确定“理想的”、能达到最优效果的ε与α的值非常困难。
技术实现要素：
4.本发明旨在解决上述问题，提供一种基于粒子分裂的sph计算优化方法及系统。
5.第一方面，本发明提供一种基于粒子分裂的sph计算优化方法，包括大粒子、小粒子和代理粒子；包括如下步骤：所述代理粒子指在优化过程中发生作用的大粒子与小粒子；所述代理粒子与在优化过程中和小粒子发生所作用的大粒子相关联；定义代理粒子是为保证粒子优化算法的计算精度，防止由于大-小粒子的直接作用产生计算误差；步骤一、将符合粒子优化准则的大粒子分裂为小粒子；所述优化准则一种为划定优化区域，一种为选取粒子类型；两种优化准则可根据实际优化需求选择一种或同时采用两种。
6.步骤二、重置代理粒子的物理参数；包括所述代理粒子的质量、粒子间距与前述大粒子分裂时产生的小粒子相同；在数值模拟过程中，当大粒子与小粒子发生作用时，大粒子的关联代理粒子将取代大粒子并与小粒子相互作用，将所述代理粒子的密度、速度、压力、光滑长度与代理粒子相关联的大粒子相同；步骤三、遍历计算域搜索粒子对，并确定粒子在当前时间步所属的计算域；
步骤四、通过求解连续性方程、动量方程及能量方程，获得大粒子、代理粒子及小粒子所承载的粒子密度dρ、粒子速度dv；步骤五、通过下式将代理粒子所承载的粒子密度dρ、粒子速度dv参量进行平均并传递给与代理粒子相关联的大粒子：式中，、分别代表过渡区大粒子及其关联的代理粒子的物理参量；n为粒子分裂的数量；一般本发明所述计算优化方法中n=ke，其中k≥2；e取值2或3，在二维情况下取值为2，三维情况下取值为3；步骤六、更新大粒子、小粒子的粒子参数，当前时间步计算结束；判断时间步次数是否大于预设次数；若否，则继续重复前述步骤进行下一时间步计算；若是，则优化结束。
7.进一步，本发明所述基于粒子分裂的sph计算优化方法，所述计算域包括三个子区域，分别为大粒子区、过渡区及小粒子区；所述大粒子区为不与小粒子发生作用的大粒子构成的区域；所述过渡区为与小粒子发生作用的大粒子构成的区域；所述小粒子区为由小粒子构成的区域。
8.进一步，本发明所述基于粒子分裂的sph计算优化方法，其特征在于：所述遍历计算域搜索粒子对时，按不同子区域搜索相应粒子对，并确定当前时间步粒子所属的子区域。粒子属于子区域的信息将用于下一时间步中的粒子对搜索及粒子信息计算。
9.进一步，本发明所述基于粒子分裂的sph计算优化方法，为保证在粒子分裂过程中的质量、动量及动能守恒，所述大粒子和小粒子的质量m、密度ρ及速度v设定为：密度ρ及速度v设定为：密度ρ及速度v设定为：；式中，mf为小粒子质量、mc为大粒子质量、ρf为小粒子密度、ρc为大粒子密度、vf为小粒子速度、vc为大粒子速度；所述大粒子与分裂形成的小粒子之间的间距、及小粒子的光滑长度hf设定为：定为：式中，dpc表示大粒子的初始粒子间距；hc为大粒子的光滑长度；α为光滑比，定义大-小粒子的光滑长度之比；ε为分离系数，定义大-小粒子间的水平和竖直距离。
10.进一步，本发明所述基于粒子分裂的sph计算优化方法，所述大粒子区内粒子的临近粒子是大粒子；大粒子区域内某一粒子i的临近粒子j位于过渡区时，j粒子以大粒子形式
参与i粒子相关物理参数的计算；所述过渡区内代理粒子的临近粒子是小粒子或大粒子代理粒子；过渡区外大粒子j位于过渡区内代理粒子i的支持域内时，j粒子的代理粒子参与计算代理粒子i的相关物理参数；所述小粒子区的小粒子为建模时指定的小粒子或由大粒子分裂产生的小粒子；小粒子的临近粒子是小粒子或过渡区的代理粒子。
11.进一步，本发明所述基于粒子分裂的sph计算优化方法，所述大粒子区内的流场信息由光滑长度为hc的核函数插值计算；所述过渡区内的流场信息通过大粒子的代理粒子计算，由光滑长度为hr=hc的核函数插值计算；所述小粒子区内的流场信息由光滑长度为hf=αhc的核函数插值计算。
12.进一步，本发明所述基于粒子分裂的sph计算优化方法，所述粒子分裂采用沿坐标轴向的正方形粒子分裂模型；所述正方形粒子分裂模型将大粒子分裂为四个小粒子，小粒子分布在以大粒子为中心的正方形的四个角上，正方形的边与坐标轴平行。
13.进一步，本发明所述基于粒子分裂的sph计算优化方法，所述正方形粒子分裂模型的分离系数ε为0.5。首先，ε=0.5有助于减小粒子分裂后小粒子之间的重叠，改善粒子秩序；其次，ε取固定值后，在数值模拟过程中，唯一需要调整的参数只有光滑比α，这大大提高了本发明所述粒子优化算法的实用性。
14.第二方面，本发明提供一种基于粒子分裂的sph计算优化系统，包分裂单元、参数重置单元、遍历搜索单元和分析输出单元；所述分裂单元用于将符合粒子优化准则的大粒子分裂为小粒子；所述参数重置单元用于重置代理粒子的物理参数；包括所述代理粒子的质量、粒子间距与前述大粒子分裂时产生的小粒子相同；所述代理粒子的密度、速度、压力、光滑长度与代理粒子相关联的大粒子相同；所述遍历搜索单元用于遍历计算域搜索粒子对，并确定粒子在当前时间步所属的计算域；所述分析输出单元用于通过求解连续性方程、动量方程及能量方程，获得大粒子、代理粒子及小粒子所承载的粒子密度dρ、粒子速度dv；并通过下式将代理粒子所承载的粒子密度dρ、粒子速度dv参量进行平均并传递给与代理粒子相关联的大粒子：；式中，、分别代表过渡区大粒子及其关联的代理粒子的物理参量；然后更新大粒子、小粒子的粒子参数，当前时间步计算结束；判断时间步次数是否大于预设次数；若否，则继续重复前述步骤进行下一时间步计算；若是，则优化结束；输出优化结果。
15.第三方面，本发明提供一种基于粒子分裂的sph计算优化系统，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其特征在于：所述处理器执行所述程序时实现第一方面所述优化方法的步骤。
16.本发明所述基于粒子分裂的sph计算优化方法及系统，根据粒子优化准则，将满足
优化准则的大粒子分裂为小粒子，将大-小粒子的接触区域定义为过渡域，引入代理粒子的概念，利用代理粒子代替大粒子与小粒子进行作用。保持了sph方法在大-小粒子的接触区域的空间插值精度，避免了由于大-小粒子的直接接触而造成的插值误差和数值震荡，在保证计算精度的同时有效提高了数值模拟效率。
附图说明
17.图1为本发明实施例所述基于粒子分裂的sph计算优化方法粒子配置示意图；图2为本发明实施例所述基于粒子分裂的sph计算优化方法计算域示意图；图3为本发明实施例所述基于粒子分裂的sph计算优化方法流程示意图；图4为本发明实施例三所述基于粒子分裂的sph计算优化系统结构示意图；图5为本发明实施例四所述基于粒子分裂的sph计算优化系统结构示意图；图6为本发明实施例二所述小球落水算例粒子配置示意图；图7为本发明实施例二所述小球入水深度的仿真与实验结果对比示意图；图8为本发明实施例二所述不同粒子配置时流场形态及压力场随时间变化过程示意图。
具体实施方式
18.下面通过附图及实施例对本发明所述基于粒子分裂的sph计算优化方法及系统进行详细说明。
19.实施例一本公开实施例采用沿坐标轴向的正方形粒子分裂模型。该粒子分裂模型将大粒子分裂为四个小粒子，小粒子分布在以大粒子为中心的正方形的四个角上，正方形的边与坐标轴平行。为保证在粒子分裂过程中的质量、动量及动能守恒，对分裂后形成小粒子的质量、密度及速度定义如下：
ꢀꢀ
（1）
ꢀꢀ
（2）
ꢀꢀ
（3）式中，mf为小粒子质量、mc为大粒子质量、ρf为小粒子密度、ρc为大粒子密度、vf为小粒子速度、vc为大粒子速度；大粒子与分裂形成的小粒子之间的间距、及小粒子的光滑长度hf设定为：
ꢀꢀ
（4）
ꢀꢀ
（5）式中，dpc表示大粒子的初始粒子间距，hc为大粒子的光滑长度，α、ε分别称为光滑比和分离系数，分别定义大-小粒子的光滑长度之比和大-小粒子间的水平和竖直距离。
20.在本公开实施例中粒子分裂模型的分离系数ε固定为0.5，这样做有以下两方面的优势：首先，ε=0.5有助于减小粒子分裂后小粒子之间的重叠，改善粒子秩序；其次，ε取固定
值后，在数值模拟过程中，唯一需要调整的参数只有光滑比α，这大大提高本实施例所述粒子优化算法实用性及效率。
21.当粒子优化算法被引入sph计算框架后，sph方法的连续性会下降。因此，为保证粒子优化算法的计算精度，防止由于大-小粒子的直接作用产生计算误差，本发明所述基于粒子分裂的sph优化算法中引入了一种新型粒子“代理粒子”，并基于代理粒子提出了一种大-小粒子的混合作用格式。
22.在本公开实施例所述基于粒子分裂的sph计算优化方法中共包括三种类型的粒子：大粒子、小粒子及代理粒子。三种类型的粒子配置如图1所示，代理粒子是与大粒子相关联的，在本公开实施例中每个大粒子都关联着四个代理粒子。因此，在本公开实施例中的代理粒子的质量及粒子间距与粒子分裂时产生的小粒子相同，即：
ꢀꢀ
（6）
ꢀꢀ
（7）式中，下标r表示代理粒子的相关物理量，mr是代理粒子的质量，dpr、dpf分别为代理粒子与小粒子的初始粒子间距。
23.在数值模拟计算过程中，当大粒子与小粒子发生作用时，大粒子的关联代理粒子将取代大粒子并与小粒子相互作用，因此规定代理粒子的密度、速度、压力、光滑长度等参数与其关联的大粒子相同，即：
ꢀꢀ
（8）
ꢀꢀ
（9）
ꢀꢀ
（10）
ꢀꢀ
（11）本公开实施例所述sph计算优化方法中，整个计算域可以划分为三个子区域，分别为：大粒子区、过渡区及小粒子区。如图2所示，（a）为大粒子区，（b）为过渡区，（c）为小粒子区，三个子区域的定义及相互作用方式如下：（a）大粒子区：由不与小粒子发生作用的大粒子构成的区域。如图2所示，该区域内粒子的临近粒子只能是大粒子。当该区域内某一粒子i的临近粒子j位于过渡区时，j粒子以大粒子形式参与i粒子相关物理量的计算。该区域内的流场信息由光滑长度为hc的核函数插值计算。
24.（b）过渡区：由与小粒子发生作用的大粒子构成的区域。如图2所示，过渡区内的流场信息通过大粒子的代理粒子计算。代理粒子的临近粒子只能是小粒子或大粒子代理粒子。当过渡区外大粒子j位于过渡区内代理粒子i的支持域内时，j粒子的代理粒子参与计算代理粒子i的相关物理量。该区域内的流场信息由光滑长度为hr=hc的核函数插值计算。
25.（c）小粒子区：由小粒子构成的区域，这些小粒子可以是建模时指定的，也可以是模拟过程中由大粒子分裂产生的。如图2所示，小粒子的临近粒子只能是小粒子或过渡区的代理粒子。该区域内的流场信息由光滑长度为hf=αhc的核函数插值计算。
26.由以上分析可知，本公开实施例所述基于粒子分裂的sph计算优化方法中，过渡区的粒子可被看作混合粒子：当与大粒子区粒子相互作用时，该区域内的粒子是大粒子；当与小粒子发生作用或过渡区内粒子自身相互作用时，过渡区内的代理粒子将代替大粒子参与流场计算；这种处理方法防止了由于大-小粒子直接相互作用而带来的计算误差和数值不稳定。因此，我们将以上不同子区域内粒子相互作用方式也称为混合粒子作用格式。
27.本公开实施例所述基于粒子分裂的sph计算优化方法如图3所示，包括如下步骤：步骤一、将符合粒子优化准则的大粒子分裂为小粒子；本实施例中采取的优化准则为划定优化区域；步骤二、重置代理粒子的物理参数；包括所述代理粒子的质量、粒子间距与前述大粒子分裂时产生的小粒子相同；在数值模拟过程中，当大粒子与小粒子发生作用时，大粒子的关联代理粒子将取代大粒子并与小粒子相互作用，将所述代理粒子的密度、速度、压力、光滑长度与代理粒子相关联的大粒子相同；步骤三、遍历计算域搜索粒子对，并确定粒子在当前时间步所属的计算域；遍历计算域搜索粒子对时，按不同子区域搜索相应粒子对，并确定当前时间步粒子所属的子区域；粒子属于子区域的信息将用于下一时间步中的粒子对搜索及粒子信息计算；步骤四、通过求解连续性方程、动量方程及能量方程，获得大粒子、代理粒子及小粒子所承载的粒子密度dρ、粒子速度dv；通过下式将代理粒子所承载的粒子密度dρ、粒子速度dv参量进行平均并传递给与代理粒子相关联的大粒子：
ꢀꢀ
（12）式中，、分别代表过渡区大粒子及其关联的代理粒子的物理参量；步骤五、更新大粒子、小粒子的粒子参数，当前时间步计算结束；步骤六、判断时间步次数是否大于预设次数；若否，则继续重复前述步骤进行下一时间步计算；若是，则优化结束。
28.本公开实施例所述sph计算优化方法从sph方法的lagrange粒子本质出发，通过对某一关键区域进行粒子优化的手段实现变分辨率数值模拟，保持了sph方法在大-小粒子的接触区域的空间插值精度，避免了由于大-小粒子的直接接触而造成的插值误差和数值震荡，在保证计算精度的同时有效提高了数值模拟效率。
29.实施例二本实施例基于前述实施例一所述基于粒子分裂的sph计算优化方法，用水球落水算例测试本发明所述计算优化方法对于包含流体与刚体相互作用及自由表面流动问题的计算效果。
30.小球落水算例的模型及粒子配置如图6所示，其中水的宽度和深度分别为2m、0.5m，小球直径r=5.5
×
10-2
m，重力加速度g=-9.81m/s2。计算开始时，小球刚接触水面，根据grennhow的实验数据，此时小球速度为v0=2.955m/s，声速cs=10v0=29.55m/s，人工粘性参数α
av
=0.05。小球落水过程的关键是小球与水之间接触力的计算，因此，本实施例中以小球为中心，设置边长为6r的粒子优化域；在计算过程中，优化域随着小球的运动而移动。水采用
40000个大粒子离散，粒子间距dpc=5.0
×
10-3
m，光滑长度hc=1.1dpc。
31.本公开实施例分别对比了以下粒子配置时的计算效果：本发明所述基于粒子分裂的sph优化算法，ε=0.5，α=0.5、0.6、0.7、0.75；reyes l
ó
pez粒子优化算法，ε=0.5，α=0.7；vacondio粒子优化算法，ε=0.4，α=0.9；传统单精度sph算法，流体粒子数n=4096、10000、40000、160000。
32.如图7所示为小球入水深度的计算结果与实验结果的对比，传统单精度sph算法的计算结果随着粒子数的增大而逐渐接近实验值。α=0.7、0.75时本发明所述粒子优化算法的计算结果与n=160000时单精度sph方法的计算结果非常接近。由于优化域随着小球的运动而移动，计算过程中的粒子数随着粒子优化而不断增加，因此，本实施例中对比了不同粒子配置情况下，计算50000个时间步的cpu总耗时（见表1）。
33.表 1 小球落水算例：不同粒子配置的计算效率对比由表1可见，本实施例与reyes l
ó
pez粒子优化算法的流体粒子数均在60000以下，与n=160000时单精度sph方法相比，本实施例与reyes l
ó
pez粒子优化算法节省了60%以上的计算时间。
34.为进一步验证粒子优化算法的计算效果，本实施例还对仿真获得的小球落水过程中的流场形态进行了分析。图8为传统单精度sph方法在流体粒子数n=160000时计算出的流场形态与压力与本发明所述方法在n=40000、α=0.7时的计算结果对比，（其中，每幅图像的左边为传统单精度sph方法使用160000个流体粒子时的仿真结果；右边为使用40000个流体粒子和新型粒子优化算法α=0.7时的仿真结果，虚线标出的是粒子优化域）可以看出，在小球落水的各个时刻，本发明所述方法与传统sph方法的计算结果在流场形态、压力波形态等方面非常一致，本发明所述方法计算出的压力场在大-小粒子区域之间过渡平滑。挡光滑比α变化时，压力场与速度场的变化一致，α=0.5的计算结果在大-小粒子接触区域存在轻微震荡，α=0.6、0.75时的计算结果与α=0.7时的计算结果基本一致。
35.由以上算例分析可见，对于小球落水问题，α=0.7时本发明所述方法的计算效率略低于reyes l
ó
pez方法，但是，α=0.6、0.7、0.75时，本发明所述方法在小球落水深度及压力场、速度场等方面的计算结果更为准确。
36.实施例三基于与前述实施例中一种基于粒子分裂的sph计算优化同样发明构思，本实施例
公开了一种基于粒子分裂的sph计算优化系统，如图4所示，所述优化系统包分裂单元、参数重置单元、遍历搜索单元和分析输出单元；其中所述分裂单元用于将符合粒子优化准则的大粒子分裂为小粒子；所述参数重置单元用于重置代理粒子的物理参数；包括所述代理粒子的质量、粒子间距与前述大粒子分裂时产生的小粒子相同；所述代理粒子的密度、速度、压力、光滑长度与代理粒子相关联的大粒子相同；所述遍历搜索单元用于遍历计算域搜索粒子对，并确定粒子在当前时间步所属的计算域；所述分析输出单元用于通过求解连续性方程、动量方程及能量方程，获得大粒子、代理粒子及小粒子所承载的粒子密度dρ、粒子速度dv；并通过下式（12）将代理粒子所承载的粒子密度dρ、粒子速度dv参量进行平均并传递给与代理粒子相关联的大粒子： （12）；式中，、分别代表过渡区大粒子及其关联的代理粒子的物理参量；然后更新大粒子、小粒子的粒子参数，当前时间步计算结束；判断时间步次数是否大于预设次数；若否，则继续重复前述步骤进行下一时间步计算；若是，则优化结束；输出优化结果。
37.通过对前述基于粒子分裂的sph计算优化方法的详细描述，本领域技术人员可以清楚的知道本实施例所述基于粒子分裂的sph计算优化系统的实施方法，为了说明书的简洁，在此不再详述。
38.实施例四基于与前述实施例中基于粒子分裂的sph计算优化方法的发明构思，本发明还提供一种基于粒子分裂的sph计算优化系统，包括处理器、存储器，所述存储器上存储有计算机程序，该程序被处理器执行时实现实施例一所述基于粒子分裂的sph计算优化方法的任意步骤。
39.图5为本实施例所述基于粒子分裂的sph计算优化系统的结构示意图。其中，总线可以包括任意数量的互联的总线和桥，总线将包括由处理器代表的一个或多个处理器和存储器代表的存储器的各种电路链接在一起。总线还可以将诸如外围设备、稳压器和功率管理电路等之类的各种其他电路链接在一起，这些都是本领域所公知的，因此，此处不再对其进行进一步描述。总线接口在总线和接收器和发送器之间提供接口。接收器和发送器可以是同一个元件，用于在传输介质上与各种其他系统通信的单元。处理器负责管理总线和通常的处理，而存储器被用于存储处理器在执行操作时所使用的数据。
40.本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质（包括但不限于磁盘存储器、cd-rom、光学存储器等）上实施的计算机程序产品的形式。
41.本发明是参照根据本发明实施例的方法、设备（系统）、和计算机程序产品的流程
图来描述的。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程中指定的功能的系统。
42.这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令系统的制造品，该指令系统实现在流程图一个流程或多个流程一个方框或多个方框中指定的功能。
43.这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程中指定的功能的步骤。
44.尽管已描述了本发明的优选实施例，但本领域内的技术人员一旦得知了基本创造性概念，则可对这些实施例进行另外的变更和修改，若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。
再多了解一些
本文用于企业家、创业者技术爱好者查询，结果仅供参考。
