[发明专利]一种基于ISPH提高飞溅流体稳定性的实现方法

摘要

本发明是一种基于ISPH(Incompressible Smoothed Particle Hydrodynamics，不可压缩光滑粒子流体动力学)提高飞溅流体稳定性的实现方法，针对传统ISPH方法在模拟扰动剧烈的流体时出现的数值不稳定问题，提出了两种改进方法，即改进泊松压力方程源项和改进拉普拉斯算子模拟。首先，提出了对泊松压力方程源项的替换方法，使得密度变化率计算更加准确；其次，提出了高阶拉普拉斯算子的构建模型，降低了离散化求解过程的数值耗散，并进一步选取高阶光滑核函数(高阶B‑样条函数)与高阶拉普拉斯算子结合使用，降低核函数二阶导数对粒子扰动的敏感性。这两种方法互相独立，却又相辅相成。实验结果表明，将这两种方法叠加使用，流体在剧烈变动时也能保证较好的压力稳定性与真实性。

主权项

1.一种基于ISPH(Incompressible Smoothed Particle Hydrodynamics，不可压缩光滑粒子流体动力学)提高飞溅流体稳定性的实现方法，其特征在于实现步骤如下：步骤一：改进泊松压力方程源项。步骤1.1，构造泊松压力方程源项求解模型。飞溅场景下的流体运动剧烈且压力变化率和变化范围都较大，所以很容易受到误差累计的影响。因此，本发明提出了改进的泊松压力方程源项求解模型，如式(1)所示。式中的*代表该表达式在预测步骤中进行计算。由于预测步骤与修正步骤都是在相同长度的时间步长Δt中进行的，而为了保证不可压缩性，预测步骤开始时候的密度ρ_start和修正步骤完成以后的密度ρ_end理论上应该都与初始密度ρ₀相等，但由于SPH方法的离散化求解方式，则三者之间存在误差相等关系，用尖括号表示，形如式(2)所示。<ρ₀>＝<ρ_start>＝<ρ_end>   (2)由于ρ_end未知，预测步骤中密度从ρ_start变化到ρ_*可以通过本发明提出的式(1)计算出，因此改进后的泊松压力源项准确性优于传统的不可压缩SPH方法的计算，在程序实验环节中也证实了这一结论。步骤1.2，推导泊松压力方程求解模型。将本发明提出的密度变化率求解公式(1)代入传统不可压缩流体的泊松压力方程中，就得到了改进后的泊松压力方程求解模型，如式(3)所示。对式(1)进行计算简化，简化过程如式(4)所示。将式(4)应用到式(3)，就得到了本发明提出的改进的泊松压力求解模型，形如式(5)所示。SPH方法的核函数为轴对称凸函数，所以根据式(5)可知，当邻近粒子趋向目标粒子移动时，密度随时间的变化率将会增大。并且目标粒子所在位置的压力也会随之增强，因而对邻近粒子施加了更大的排斥力。并且泊松压力方程中的源项的变化与预测过程中的速度增量Δu_**成正比。所以，本发明提出的改进的泊松压力方程求解模型不会影响流体的不可压缩性。步骤二：改进拉普拉斯算子模拟。步骤2.1，改进高阶核函数拉普拉斯算子。本发明在不可压缩的SPH方法的基础上引入直接二阶微分法并且采用高阶B‑样条核函数，来改进SPH方法，使得在模拟飞溅等压力变动较大的场景下帧率更加稳定。通过分析SPH方法中梯度模型的散度，容易得出如式(6)所示的拉普拉斯算子的一阶插值公式。尖括号代表该项存在误差相等关系。其中，φ表示任意物理量，w表示核函数，m表示粒子的质量。下标i和j分别对应当前粒子和当前粒子的邻近粒子。而由于不可压缩性的流体微元在运动过程中由质量守恒可知其体积V不会发生变化，因而在二维场景下可以推导出直接二阶微分方法的拉普拉斯算子表达式，形如式(7)。步骤2.2，求解泊松压力方程。利用引入的拉普拉斯算子的直接二阶微分方法来求解不可压缩SPH方法中的泊松压力方程，其过程如式(8)所示。由于直接二阶微分方法在对拉普拉斯算子进行离散求解时，并没有真正消除或者降低求导次数，而是将二阶导数计算转移给了核函数。因此，当核函数的阶数低于3阶时，二阶求导计算将对其性质造成严重影响，从而使得求解结果对粒子的扰动过于敏感。而本发明的实现方法充分考虑到了这一缺陷，选择B‑样条函数作为SPH方法的核函数，如式(9)所示。式中q＝r/h，r表示粒子间的距离，h表示光滑长度。在实际实现中，本发明提出的高阶拉普拉斯构造和高阶光滑核函数相结合的方法，在飞溅流体场景下，稳定性和真实性优于改进前的标准不可压缩SPH方法。