CN109783935A - 一种基于isph提高飞溅流体稳定性的实现方法 - Google Patents

Abstract

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

Description

一种基于ISPH提高飞溅流体稳定性的实现方法

技术领域

本发明涉及计算机图形学领域，涉及一种基于ISPH(Incompressible SmoothedParticle Hydrodynamics，不可压缩光滑粒子流体动力学)提高飞溅流体稳定性的实现方法，具体包含两种方法，即改进泊松压力方程源项和改进拉普拉斯算子模拟。

背景技术

流体的飞溅现象是自然界中最常见的现象之一，常见于瀑布、浪花拍岸等场景中。在计算机图形学中，模拟流体的飞溅现象最重要的就是模拟出该现象场景的真实感。有两种基本的思路来增强模拟飞溅流体时的真实性，一种是人工加入相关场景的方法，比如人工加入泡沫和气泡贴图，让飞溅看起来更真实；另一种则是通过提高SPH(SmoothedParticle Hydrodynamics，光滑粒子流体动力学)方法模拟流体时的稳定性和真实性，从而可以从流体运动的真实性自然延伸至飞溅现象真实性。

SPH方法是一种用于求解连续介质动力学方程的无网格方法，它通过引入空间场函数和核函数的概念，将流体控制方程离散化。该方法避免了网格方法中存在的网格缠绕和扭曲问题，在计算机物理动画中主要用来进行流体的运动模拟。

由于飞溅场景中流体的密度和压力变化率都非常大，因此该场景对于离散化求解的准确性要求较高。传统的不可压缩SPH方法尽管在很多场景下都能有非常高效的模拟效果，例如电影特效等，但是由于数值耗散的存在，一般不能长时间准确模拟流体的运动变化，在对数值准确性要求较高的场景下还有很大的改进空间。

针对传统方法的不足之处，本发明将给出整个模拟过程中压力求解的流程，提出大幅度提高计算准确性的新方法，使得对飞溅性流体的模拟效果更加稳定与真实。

在使用SPH方法模拟飞溅性不可压缩流体方面，国内外很多大学和研究机构都进行了全面而深入的研究，并取得了突破性进展。

1996年，S.Koshizuka等提出了移动粒子半隐式方法(Moving Particle Semi-implicit method，MPS)来改进SPH方法，以逼真模拟不可压缩流体。在MPS方法中，用于粒子相互作用模型的权重函数得到了改进，从而提高了MPS方法的数值稳定性。但是，MPS方法需要额外的计算时间来生成邻居列表，其中记录了相邻的粒子数量。每次生成列表都要遍历计算所有粒子，时间复杂度为I(N²)，所以当进行大规模模拟时，时间开销非常巨大(参考文件1：S.Koshizuka,Y.Oka.Moving particle semi-implicit method for fragmentationof incompressible fluid[J].Nuclear Science and Engineering,1996,123(3):421-434)。

1998年，S.Koshizuka等提出对MPS方法的改进策略。其中提出一种算法来将时间复杂度减少到N^1.5。并通过使用网格的数值模拟方法来进行自由表面的模拟。但自由表面经常发生碎裂和折叠现象，使用传统网格方法需要存储大量顶点及网格的位移数据，无法满足实时的帧率要求。其他研究文献中提到用静止网格分析断裂波，将自由表面用多段组成的多边形来表示，在每个时间步长中计算多边形的运动。(参考文件2：KoshizukaS,OkamotoK,FurutaK.Development of computational techniques for nuclearengineering[J].ProgressinNuclearEnergy,1998,32(1-2):209-222)

2003年，Songdong Shao等提出了一个严格不可压缩的SPH模型来模拟自由表面流动。其使用基本方法类似于Cummins等在SPH预测方法中提到的对没有自由表面的流体的模拟方式。它采用严格不可压缩的SPH约束方程，其中压力不是显式通过求解状态方程获得的，而是通过求解压力泊松方程来获得的。这种对SPH方法的扩展也进一步说明了SPH方法的灵活性。(参考文件3：Shao S,Lo E Y M.Incompressible SPH method for simulatingNewtonian and non-Newtonian flows with a free surface[J].Advances in WaterResources,2003,26(7):787-800.)

以上方法提供了对流体表面的模拟基础。与此同时，近年来，更多的研究人员开始重视如何将真实感的流体表面延申到高飞溅流体细节模拟。

2007年，N.Thürey等提出了在流体中加入气泡和泡沫的方法。这种方法是通过浅水方程与粒子法模拟的气泡和泡沫耦合来实现的。并通过设计球形漩涡，从而可以在气泡周围产生具有漩涡的气流场。(参考文件4：Thürey N,Sadlo F,Schirm S,et al.Real-timesimulations of bubbles and foam within a shallow water framework[J].2007.)

2017年，Jong-Hyun Kim等提出了一种基于物理的有效方法来模拟飞溅效果。其基本原理是在流体周围设计小的空气腔，并在计算表面粒子密度时采用在圆形区域内统计粒子个数与种类。在模拟瀑布与海浪等飞溅很大的场景时，这种方法取得了非常逼真的效果。(参考文件5：Kim J H,Kim W,Lee J.Physics-inspired approach to realistic andstable water spray with narrowband air particles[J].Visual Computer,2017:1-11.)

Jena等使用MPS无网格粒子法研究了液体在容器中的晃动行为，通过正弦波和震动对内部液体进行激励以形成飞溅等不稳定行为，然后检测数值模型的有效性。在该论文的实验中，与高频率震动相比，低频震动可以产生更高的动态响应幅度，如冲击压力载荷和飞溅对容器上壁的作用。(参考文件6：Jena D,Biswal K C.Violent Sloshing and WaveImpact in a Seismically Excited Liquid-Filled Tank:Meshfree Particle Approach[J].Journal of Engineering Mechanics,2018,144(3):04017182.)

以上方法在流体平缓运动的场景下效率极高；但在流体飞溅运动的场景下，这些方法会造成压力稳定性不够和场景失真。这是由于，传统方法在构建泊松压力方程的过程中，密度变化率离散化不精确，从而产生累计误差导致的。

本发明提出一种飞溅性流体建模的方法，通过改进泊松压力方程源项与拉普拉斯算子的离散化方法，从而更加精确地求解飞溅性流体物理模型。

这里先对基础理论加以说明。

式(1)描述了不可压缩流体的质量守恒，最初是在欧拉方程及求解法中提出的。可压缩流体的质量守恒方程也被称为连续性方程，在密度场中求解物质导数便可得到该方程的一般形式，如式(1)所示。

其中ρ表示流体的粒子数密度，t表示时间，u表示流体的速度。运算符·为空间向量的点乘，表示空间向量对空间维度的偏导，即又称为梯度。

式(1)左侧第一项为粒子质量随时间的变化率，第二项为例子边界上流入边界内的粒子质量。由于这里描述的流体是不可压缩的，因而质量密度在流体运动时保持不变，即有和

下面依次介绍传统ISPH方法的求解步骤。

二维场景下的自由表面流体流动的控制方程由质量守恒方程和动量守恒方程组成。本发明研究粒子方法的流体模拟，因此将控制方程表示为拉格朗日形式，质量守恒方程和动量守恒方程分别如式(2)和式(3)所示。式(2)是由式(1)两边同时乘以得到。式(3)则由纳维-斯托克斯(Navier-Stokes，N-S)方程得来。

其中ρ表示流体的粒子数密度，t表示时间，u表示流体的速度，P表示压力，g表示重力加速度，τ表示剪切应力(包括黏性力)。由于保留了黏性力特征，所以式(2)和式(3)既适用于牛顿型流体(粘度系数为常数，如空气、水等)也适用于非牛顿型流体(粘度系数不是常数，如聚合溶液、含有悬浮杂质或纤维的流体等)。

式(2)描述了可压缩流体的质量守恒形式。在不可压缩流体的计算中，引入可压缩流体的质量守恒公式来预测出每个时间步长Δt中流体密度的偏差，然后在修正步骤中使用密度偏差值来强制保证流体的不可压缩性。也就是说，ISPH方法在传统SPH方法的基础上，加入了预测步骤和压力修正步骤，之后根据得到的泊松方程源项，求解泊松方程，得到粒子的属性。因此，上述两个关键步骤影响了最终不可压缩水体的流动特性。下面对这两个步骤进行介绍。

1、速度与密度预测

流体粒子的内压P在SPH方法的计算中由粒子的内密度所决定。因此内压与内密度之间有直接的相互作用关系，修正密度将会直接影响内压力的值。在Navier-Stokes中如果只考虑应力和引力，则可以由式(4)、式(5)和式6)得到预测与修正步骤的中间时刻粒子的速度和位置。

u_*＝u_t+Δu_* (5)

r_*＝r_t+u_*Δt (6)

这里u_t和r_t分别表示粒子在时刻t的速度和位置。u_*和r_*表示预测步骤结束后的粒子的中间速度和位置。Δu_*表示预测过程中粒子速度的改变量，Δt表示时间步长值。式(4)是由式(3)移项变形得到的。式(5)描述了时间步长Δt内的中间速度与当前时刻速度的关系。式(6)描述了步长时间内的粒子速度与位置的关系。

2、压力-密度修正

ISPH是通过粒子密度来校正其速度和压力属性的。第一个步骤的计算并不能保证流体的不可压缩性。r_*位置处的流体密度ρ_*在SPH方法中由式(7)给出。

ρ_*＝∑_jm_jW(|r_t-r_j|,h) (7)

式中r_t表示预测步骤开始之前的粒子位置坐标，r_j和m_j分别表示目标粒子的邻近粒子和邻近粒子的质量。从式(7)可以看出，ρ_*与时刻t的粒子质量、光滑函数相关，并不是一个恒定的值，它与不可压缩流体密度ρ₀之间存在偏差。

因此，修正步骤要做的操作就是把流体密度调整为初始值。之后利用密度对压力进行校正。在计算过程中，压力项将被用于更新当前步长中间步骤中计算出的粒子速度。至此完成对速度的校正。相关步骤如式(8)和式(9)。

u_t+1＝u_*+Δu_** (9)

这里的Δu_**表示校正过程中粒子速度的改变量，ρ_*表示预测步骤结束后与修正步骤开始前密度的中间量，P_t+1和u_t+1表示在t+1时刻的流体压力和粒子速度。

最后，选取足够小的时间步长Δt，使用平均速度来近似计算出粒子的位置，如式(10)所示。

其中r_t和r_t+1分别表示t和t+1时刻粒子的位置。

执行强制不可压缩性所需的压力项来自于式(2)，即质量守恒方程。形如式(11)所示。

这里的ρ₀表示初始密度值，也是预测修正过程完成后的最终密度值。至此，传统的ISPH预测与修正步骤结束。

结合式(8)和式(11)，就得到了压力的泊松方程，形如式(12)所示。其中就是柏松压力方程源项。源项是一个广义量，它代表了那些不能包括到控制方程的非稳态项、对流项与扩散项中的所有其它各项之和。通过求解该方程，最终更新粒子位置属性，完成当前步长的迭代。

综上所述，ISPH方法通过预测与修正步骤来保证不可压缩性，其核心就是通过假设密度还原来求解出粒子内压然后隐式求解出速度修正量Δu_**。因此，如何精确地离散化求解出泊松压力方程，将直接影响整个修正过程的精确度，从而影响流体的不可压缩性。

基于上述背景，本发明对泊松方程源项以及拉普拉斯算子(即二阶微分项，形如)进行改进，从而大幅提高飞溅流体的稳定性与真实性。

式(12)中，修正时间内密度的变化率就是泊松方程源项。基于SPH方法，使用微分插值来将连续的复杂方程离散线性化。在模拟不可压缩流体时，对粒子密度ρ的修正是整个算法的核心。粒子i的密度ρ_i由其邻近粒子质量累加和得到，如式(13)所示。

ρ_i＝∑_jm_jW(|r_i-r_j|,h)＝∑_jm_jW_ij (13)

式中下标i和j分别表示目标粒子i和目标粒子的邻近粒子j。m表示质量，r表示位置坐标，W表示核函数，h表示光滑长度。ρ是不特指i粒子或t时刻的粒子密度。ρ₀是步长初始时刻粒子的密度。ρ_*是为了求解步长内的粒子密度，而设置的粒子的临时密度。从ρ₀到ρ_*变化的过程是非线性的。

求得密度ρ之后，就可以通过密度求得下一时间步长上的压力p，引入反对称离散化形式来表示，形如式(14)。

采用反对称形式求解压力梯度，可以保证目标粒子i与邻近粒子j之间计算出的值跟目标粒子j与邻近粒子i之间计算出来的值方向相反、大小相等。用反对称表示法来表示力和加速度具有非常强的物理意义，可以保证粒子系统满足牛顿第三定律和动量守恒定律。

式(14)描述了传统SPH方法中压力的计算方法。而在对不可压缩流体的求解过程中，需要通过求解泊松方程来得到压力，并在这个过程中隐式求解出速度修正量Δu_**。在修正步骤中假定流体的不可压缩性都可以被完全修正。也就是说，每个中间步骤所计算出的临时密度ρ_*都可以被精确地修正为初始密度ρ₀。由于Δt足够小，且流体密度在理论上具有连续变化性，即ρ对时间t可微，则可通过线性变化代替修正时间内密度的变化率在预测步骤中构造为因此式(11)可表示为式(15)的形式。

但是，由于上述所有方程均采用粒子系统离散化来求解，所以很难保证所有目标粒子在每个时间步长结束时的密度都被修正成了初始密度ρ₀。使用ρ₀来计算密度随时间的变化率Dρ/Dt将会产生累计误差，最终造成流体模拟的压力不稳定性。因此，密度随时间的变化率Dρ/Dt的构造方法不应当依赖于ρ₀。

流体不可压缩性的满足不仅取决于泊松压力方程源项的精度，而且还取决于时间上所用数值格式的精度、积分以及微分算子的离散化，如拉普拉斯算子的离散化。拉普拉斯算子就是二阶微分算子。

接下来，对拉普拉斯算子及其三种常见离散化方法进行简单介绍。

拉普拉斯算子是水体仿真建模研究领域最常用到的运算操作之一，其在欧几里得空间中定义为二阶微分算子。对函数f运用拉普拉斯算子表示先对f做梯度运算再对梯度做散度运算例如，在流体动力学理论中，粘滞力与流体粒子速度的二阶微分成正比。SPH方法中求解粘滞力时就需要用到拉普拉斯算子，并且一般用典型的SPH微分离散法对拉普拉斯算子进行求解，如式(16)所示。

式中表示粘滞力，μ表示黏度系数。

在SPH方法研究中需要构造拉普拉斯算子时，有三种不同的方法来将连续的数学运算函数离散化表示：双求和方法、直接二阶微分方法和混合方法。其中双求和方法的基本原理是将一阶导数的差分近似代入到外层的一阶导数差分近似中。因此其包含了两次基于SPH方法的导数近似，求解过程形如式(17)所示。尖括号代表该项存在误差相等关系。

式(17)表示任意物理量φ对位置坐标x采用双求和方法近似拉普拉斯运算。

直接二阶微分法建立在对SPH插值方法的两次微分上，因此本质上是对核函数直接进行二阶微分运算，其求解过程形如式(18)所示。尖括号代表该项存在误差相等关系。

而混合方法则是兼具了前两种方法的思想，利用SPH有限差分来构建外层一阶导数，而对内层一阶导数则直接求导。双求和方法与混合方法都避免了直接对核函数进行复杂的二阶微分计算，但是在扰动较大的物理场景下，可能会扩大压力求解的误差，造成模拟失真的现象。

传统的不可压缩SPH方法在处理泊松压力方程的拉普拉斯算子时，一般使用混合方法来构建。因此拉普拉斯算子被分解成了两部分，标准SPH方法求解一阶导数和有限差分近似法求解另一个一阶导数。如式(19)所示。

式中P_ij＝P_i-P_j、r_ij＝r_i-r_j分别为顶点i和j的压强差和位置差。

发明内容

本发明提出一种飞溅性流体建模的方法，通过改进泊松压力方程源项与拉普拉斯算子的离散化方法，从而更加精确地求解飞溅性流体物理模型。

具体来讲，首先基于SPH方法，采用微分插值法求解校正步骤中的密度改变量，采用直接对核函数求导法求解密度变化率，提高泊松方程源项(源项指一个广义量，它代表了那些不能包括到控制方程的非稳态项，对流项与扩散项中的所有其它各项之和)的精确度。

此外，在简单基于SPH方法显式求解泊松压力方程中速度校正项时，由于对拉普拉斯运算采用参数近似和二阶离散方法求解，易引起数值插值计算的不稳定，且其参数调整复杂，初始化效率低。针对此不足，本发明提出一种高阶拉普拉斯离散化方法，将二阶导数运算全部转移到核函数运算中，并进一步选择高阶核函数，从而提高计算过程中的稳定性和准确性。

本发明具体包括如下步骤：

步骤一：改进泊松压力方程源项。

步骤1.1，构造泊松压力方程源项求解模型。

飞溅场景下的流体运动剧烈且压力变化率和变化范围都较大，所以很容易受到误差累计的影响。因此，本发明提出了改进的泊松压力方程源项求解模型，如式(20)所示。

式中的*代表该表达式在预测步骤中进行计算。由于预测步骤与修正步骤都是在相同长度的时间步长Δt中进行的，而为了保证不可压缩性，预测步骤开始时候的密度ρ_start和修正步骤完成以后的密度ρ_end理论上应该都与初始密度ρ₀相等，但由于SPH方法的离散化求解方式，则三者之间存在误差相等关系，用尖括号表示，形如式(21)所示。

<ρ₀>＝<ρ_start>＝<ρ_end> (21)

由于ρ_end未知，预测步骤中密度从ρ_start变化到ρ_*可以通过本发明提出的式(21)计算出，因此改进后的泊松压力源项准确性优于传统的不可压缩SPH方法的计算，在程序实验环节中也证实了这一结论。

步骤1.2，推导泊松压力方程求解模型。

将本发明提出的密度变化率求解公式(20)代入传统不可压缩流体的泊松压力式(12)中，就得到了改进后的泊松压力方程求解模型，如式(22)所示。

对式(20)进行计算简化，简化过程如式(23)所示。

将式(23)应用到式(22)，就得到了本发明提出的改进的泊松压力求解模型，形如式(24)所示。

SPH方法的核函数为轴对称凸函数，所以根据式(24)可知，当邻近粒子趋向目标粒子移动时，密度随时间的变化率将会增大。并且目标粒子所在位置的压力也会随之增强，因而对邻近粒子施加了更大的排斥力。并且泊松压力方程中的源项的变化与预测过程中的速度增量Δu_**成正比。所以，本发明提出的改进的泊松压力方程求解模型不会影响流体的不可压缩性。步骤二：改进拉普拉斯算子模拟。

步骤2.1，改进高阶核函数拉普拉斯算子。

本发明在不可压缩的SPH方法的基础上引入直接二阶微分法并且采用高阶B-样条核函数，来改进SPH方法，使得在模拟飞溅等压力变动较大的场景下帧率更加稳定。

通过分析SPH方法中梯度模型的散度，容易得出如式(25)所示的拉普拉斯算子的一阶插值公式。尖括号代表该项存在误差相等关系。

其中，φ表示任意物理量，w表示核函数，m表示粒子的质量。下标i和j分别对应当前粒子和当前粒子的邻近粒子。而由于不可压缩性的流体微元在运动过程中由质量守恒可知其体积V不会发生变化，因而在二维场景下可以推导出直接二阶微分方法的拉普拉斯算子表达式，形如式(26)。

步骤2.2，求解泊松压力方程。

利用引入的拉普拉斯算子的直接二阶微分方法来求解不可压缩SPH方法中的泊松压力方程，其过程如式(27)所示。

由于直接二阶微分方法在对拉普拉斯算子进行离散求解时，并没有真正消除或者降低求导次数，而是将二阶导数计算转移给了核函数。因此，当核函数的阶数低于3阶时，二阶求导计算将对其性质造成严重影响，从而使得求解结果对粒子的扰动过于敏感。而本发明的实现方法充分考虑到了这一缺陷，选择B-样条函数作为SPH方法的核函数，如式(28)所示。

式中q＝r/h，r表示粒子间的距离，h表示光滑长度。

在实际实现中，本发明提出的高阶拉普拉斯构造和高阶光滑核函数相结合的方法，在飞溅流体场景下，稳定性和真实性优于改进前的标准不可压缩SPH方法。

附图说明

图1是基于ISPH提高飞溅流体稳定性的实现方法流程图。

图2是飞溅性流体实验系统模块图。

图3是自由表面形变压力稳定性和准确性测试效果图。

图4是飞溅水花处压力稳定性和准确性测试效果图。

图5是三维场景构建与内部压力变化示意图。

图6是三维场景边界压力变化示意图。

图7是三维场景飞溅水花压力分布及效果图。

图8是飞溅水体准确还原稳定状态示意图。

图9是基于CPU的传统算法与本发明算法帧率对比图。

图10是基于CUDA的传统算法与本发明算法帧率对比图。

具体实施方式

下面将结合附图和实施例对本发明的技术方法及结果作进一步的详细说明。

本发明提出一种改进的泊松压力方程求解方法，求解基于ISPH的飞溅性流体建模方法。基于ISPH方法，采用微分插值法求解校正步骤中的密度改变量，采用直接对核函数求导法求解密度变化率，提高泊松方程源项的精确度。

同时，本发明提出一种高阶拉普拉斯离散化方法，将二阶导数运算全部转移到核函数运算中，选择高阶核函数，提高计算过程中的稳定性和准确性。如图1所示，下面具体对实现步骤进行说明。

步骤一：不考虑压力作用，求得速度预测值Δu_*。

步骤二：遍历邻近粒子，求解当前粒子的密度。

步骤三：遍历邻近粒子，求解核函数梯度累加和，从而计算出密度变化率。

步骤四：求解改进后的泊松压力方程，得到速度校正量Δu_**。

步骤五：通过速度变化来求解出粒子位置变化。

图2展示了飞溅性流体实验系统结构框图。实验系统主要的功能有初始化配置、邻近粒子搜索、SPH光滑核函数计算、泊松压力方程的求解、高阶拉普拉斯算子离散化应用、碰撞检测与处理、追踪表面粒子、渲染流体表面。其中初始化配置的作用是配置构建流体运动时所涉及的常数参数，以及流体粒子信息的初始化。根据实验系统所需完成的主要算法功能，可以将实验系统划分为初始配置模块、SPH处理模块、流体表面追踪模块和流体表面渲染模块。

其中初始配置模块主要实现流体模拟信息初始化功能；SPH处理模块主要完成流体模拟时粒子信息更新的计算：粒子位置的更新计算(泊松压力方程求解、碰撞检测与处理)、粒子速度的更新(强制不可压缩)、粒子存储温度更新等；流体表面追踪模块主要是完成流体表面粒子标记计算；流体表面渲染模块主要是基于Marching Cubes算法借助OpenGL工具，对流体表面进行渲染。

本发明提出了两种对泊松压力方程的解决方法，分别是高阶拉普拉斯构造方法和泊松方程源项替换方法。两者是独立的解决方法，两者与传统方法相比效率都略有降低，但都可以提高数值精度。两种方法相辅相成，通过改进泊松方程压力项使得流体有更稳定的压力分布；通过改进拉普拉斯算子，可以显著减少噪点。这两种方法都是为了提高流体压力波动较大时其变化率的稳定性，从而更加准确地模拟出飞溅等场景。图3验证了在流体自由表面发生形变时，本发明所提出的两种解决方法的有效性。

图3中(a)图显示了在流体模拟系统界面中不勾选任何改进方法，直接模拟二维场景时，自由表面形状及压力变化情况。观察可知，(a1)、(a2)中均出现了压力噪点(小白点)，这是由于泊松方程求解误差所引起的，尽管在模拟相对平静的流体时不会有太大影响，但是如果流体长时间处于扰动较大的状态下，这种误差会逐渐积累使得模拟失真。

(b)图应用了本发明提出的泊松压力方程源项改进方法，可以看到压力分层稳定性与流体表面稳定性都有所提高，没有出现白色噪点。

(c)图应用了本发明提出的高阶拉普拉斯构造法和高阶核函数相结合的方法，尽管对于压力的具体计算过程没有改进，但是提升了整个粒子系统的离散化精确度，压力噪点得到控制，但是右下角处的压力依然有压力失真现象。

(d)图中，同时使用本发明提出的两种改进方法。结果表明，形变处压力与流体下方压力模拟的真实性都有所提高。

为了验证大飞溅时本实验程序的稳定性与有效性，在系统界面中选择外界扰动大的选项，则可以得到如图4所示的场景。在程序中增大了初始外力并降低了容器上壁的高度，使流体能够受到各个方向上的撞击。

图4中(a1)与(a2)描绘了流体在撞击容器右上角时瞬时压力的分布情况。其中(a1)为本发明提出的方法，(a2)为SPH传统方法。分别对比两张图中的左上角和右上角可知，应用本发明提出的飞溅模拟方法后，压力噪点明显减少。(b1)与(b2)描绘了流体形成卷曲水花时压力的分布情况，从两张图右侧可以看出，(b1)图中流体与容器壁撞击处的压力分布明显比(b2)图中要更加真实。因此可以直观地验证出，本发明提出的基于改进泊松压力源项与高阶离散化拉普拉斯算子求解模型的算法，在飞溅场景下能够取得更加真实且更加稳定的压力分布状态。

流体运动的实验效果如图5～图7所示。

图5中第一幅图为在没有碰撞到容器壁时压力的分布，第二幅图为粒子状态下水体，第三幅图为渲染后水体效果。图6中第一幅图为水体碰撞容器壁造成的压力变化，第二幅图为粒子状态下水体，第三幅图为真实渲染后水体效果。图7中第一幅图为流体发散飞溅时压力的分布，第二幅图为粒子状态下水体，第三幅图为渲染后水体效果。从图5～图7中可知，渲染后的水体符合自然界的运动规律，并且本程序还能保证长时间运行时的稳定性与准确性。这充分说明了本发明提出的流体模拟算法的正确性。

传统的不可压缩SPH方法，在很多场景下都能够高效地模拟出流体的运动，比如电影特效和虚拟融合场景。但是在需要精确模拟时传统方法还有很大的改进空间，这也体现了SPH方法的可扩展性。例如大坝泄洪时水流对侧壁的冲击压力的计算与高飞溅水体的真实模拟，都需要更加真实的模拟算法。

图8描绘了本发明中三维仿真程序模拟出的在两个水体状态下的变化过程，左图为水体碰触容器壁回弹的渲染后效果，右图为水体平静后的渲染效果。可以看出，本发明模拟的飞溅流体非常真实。

下面，通过实验数据对比改进后的求解模型与传统不可压缩SPH方法的模拟帧率。对于含有11520个粒子的三维海浪模拟实验(图5～图7)，分别在基于CPU和基于NVIDIACUDA的平台上进行了9次独立实验，并记录了实验结果，如表1和表2所示。表1为基于CPU的实验测试数据表，表2为基于CUDA的实验测试数据表。

表1基于CPU的实验测试数据

图9是根据表1中的数据绘制得到的。从图中可以看出，传统的方法与本发明改进的方法本发明的测试机环境下，帧率都低于25fps，因此都不能实现实时模拟的效果。

表2基于CUDA的实验测试数据

图10是基于CUDA平台的实验数据绘制出的，可以看出两种算法都可以保持在300fps以上，超过了25fps的实时性标准。因此在NVIDIA CUDA GeForce 1060 6GB上运行时，两者模拟的实时性没有本质差别。

综上所述，本发明提出的基于ISPH提高飞溅流体稳定性的实现方法虽然在各计算平台上帧率略低于传统算法，但是模拟效果却有显著提升。实验结果表明，该方法不但提高了模拟的稳定和真实性，且在非极端情况下没有对程序模拟的实时性造成影响。

Claims (3)

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方法。

2.根据权利要求1所述的一种基于ISPH提高飞溅流体稳定性的实现方法，其特征在于，步骤一所构建的泊松压力方程源项，方法是：重新构造泊松方程压力源项，从而使得密度变化率计算更加准确。具体是：首先计算除压力外的其他力，并计算粒子的预测密度与预测速度；然后将密度输入到构造泊松压力方程源项中进行求解，得到速度校正量；最后通过速度变化来求解出粒子位置变化。

3.根据权利要求1所述的一种基于ISPH提高飞溅流体稳定性的实现方法，其特征在于，步骤二所述的改进拉普拉斯算子模拟，具体是：提出高阶拉普拉斯算子的构建模型，降低离散化求解过程的数值耗散；并进一步选取高阶光滑核函数(高阶B-样条函数)与高阶拉普拉斯算子结合使用，降低核函数二阶导数对粒子扰动的敏感性。