[发明专利]一种基于IISPH提高不可压缩水体模拟效率的实现方法

摘要

本发明是一种基于IISPH(Implicit Incompressibility Smooth particle hydrodynamics，隐式不可压缩光滑粒子流体动力学)方法提高不可压缩水体模拟效率的实现方法。首先将速度无散度模型引入IISPH，修正IISPH方法每个时间步长内速度场的散度误差；其次，通过改进的密度恒定模型，使得密度恒定模型与速度无散度模型共享计算因子，实现不可压缩水体的高效模拟。通过引入无散度速度模型，实现了散度误差的修正，解决了密度误差随时间不断增加的问题；通过改进密度恒定模型，实现了速度无散度模型与密度恒定模型的计算因子共享，解决了冗余计算的问题。本发明能够达到17k粒子规模下30.10帧/秒的实时模拟速度，与PCISPH(Predictive‑Corrective Incompressible SPH，基于预测‑校正的SPH)方法相比，模拟效率提高了24.17％；与IISPH方法相比，模拟效率提高了16.21％。

主权项

1.一种基于IISPH(Implicit Incompressibility Smooth particle hydrodynamics，隐式不可压缩光滑粒子流体动力学)提高不可压缩水体模拟效率的实现方法，其特征在于实现步骤如下：步骤一：将速度无散度模型引入IISPH，构建不可压缩水体模拟框架首先，基于IISPH搭建不可压缩水体模拟框架。IISPH分为预测与校正两个阶段，预测阶段使用除压力外的其他力计算粒子密度以及速度，校正阶段通过施加压力实现水体的不可压缩性。其次，改进IISPH的密度恒定模型，使其与速度无散度模型共享计算因子，减少冗余计算，提高水体模拟效率。最后，将速度无散度模型引入不可压缩水体模拟框架中。针对IISPH求解出的速度存在散度误差的问题，将速度无散度模型引入IISPH方法中，修正粒子的散度误差。步骤二：引入速度无散度模型，以修正IISPH方法的散度误差IISPH修正的粒子速度存在散度误差，即不满足本发明通过引入速度无散度模型，并使用该模型计算压力项，从而实现对粒子速度的散度误差进行修正，粒子i所受压力的大小由式(1)决定。其中表示粒子i所受的压力，表示水体粒子的体积，表示压力梯度。引入表示对粒子i的速度v进行修正的刚性系数，可得压力梯度的计算式如式(2)所示。其中m_j表示邻近粒子j的质量，表示核函数的梯度，即为满足速度无散度条件的刚性系数。粒子在运动的过程中满足动量守恒。对粒子i进行受力分析可知，粒子i所受的压力为邻近粒子j的压力之和，并且粒子间的压力是一对相互作用力。由于相互作用力大小相等，方向相反，因此，粒子i所受的压力与粒子i对邻近粒子j的压力满足其中表示粒子i对粒子j的压力，这意味着压力作为一种粒子间的作用力，以内力的形式存在，并且所有粒子的压力之和为0。从水体的整体看来，该压力不产生水体本身动能的改变，满足动能守恒定律。粒子i对邻近粒子压力的大小与粒子间的距离有关，其求解公式如式(3)所示。其中x_j表示粒子j的三维坐标，求解器必须通过改变压力值修正粒子的速度，从而使得粒子的速度满足无散度性。密度相对于时间的变化率如式(4)所示。其中W_ij表示光滑核函数，要保证水体的不可压缩性，即密度不发生变化，则式(4)求得的密度相对于时间的变化率应为0。式(4)中v_i‑v_j表示粒子间的相对速度差，该速度差由压强差导致，不同粒子受到的压力存在差别，因此导致粒子间的速度不一样。由压力引起的速度变化量如式(5)所示。其中Δv_i表示由压力引起的速度变化量，将式(5)代入式(4)，得到式(6)。由式(1)可以计算出粒子所受压力的大小，结合式(4)密度变化率为0，即可求解出刚性系数；使用该刚性系数修正粒子速度，能够使得速度的散度为0，从而使得密度相对于时间的变化率为0。将式(2)代入式(1)，得再代入式(6)，得到式(7)。根据式(7)求解出刚性系数如式(8)所示。其中且α_i的大小只与粒子位置有关。式(8)为速度无散度模型，该模型通过求解刚性系数并使用该刚性系数计算压力，从而保证水体在运动的过程中密度变化率不发生波动，即保证水体的不可压缩性。最后根据式(9)计算粒子受到的合力，并使用该合力重新计算粒子的速度。其中表示粒子i所受合力，表示除压力外其他力的合力，表示粒子i受到的邻近粒子j的压力。步骤三，改进密度恒定模型，使其与速度无散度模型共享计算因子，以修正IISPH的密度误差。在速度无散度模型中，使用满足速度无散度性的刚性系数修正粒子速度，从而保证密度变化率小于给定的散度误差值η^div。在进行无散度速度修正时需要计算α因子，本发明提出的密度恒定模型复用α计算因子，对水体粒子的密度进行进一步的修正。本发明通过密度变化率以及时间步长计算密度变化量，则t+Δt时刻的密度值ρ′_i如式(10)所示。将速度无散度模型的密度变化率求解公式(6)代入式(10)，并结合密度恒定条件(ρ_i＝ρ₀)，得到Δt时间步长内的密度误差ρ′_i‑ρ₀，如式(11)所示。其中ρ₀表示粒子初始密度，将式(1)、(2)代入式(11)。由式(12)得到刚性系数如式(13)所示。其中为满足水体密度恒定的刚性系数，α为速度无散度模型中的计算因子，在同一个时间步长内，该计算因子仅计算一次，故可以减少重复计算导致的时间开销。式(15)为密度恒定模型，该模型通过求解刚性系数并使用该刚性系数计算出的压力，从而保证水体在运动的过程中密度变化量不发生波动，即保证水体的不可压缩性。由式(11)可知，密度恒定模型修正过后粒子所受的合力如式(14)所示。使用矫正过后的粒子所受的合力对粒子的速度进行修正，结合牛顿第二定律可得修正后t+Δt时刻粒子速度v′_i的求解公式，如(15)所示。根据上述计算得出的以及式(15)对粒子的速度v_i进行修正，得到修正后的粒子速度v′_i。步骤四：利用步骤二、三所求解的物理模型，构建不可压缩水体运动场景，输出结果并显示。步骤二、三中求解的不可压缩水体模型，仅仅是基于数学物理求解得出的结果，从视觉上描述，是由一系列粒子构成，并不能逼真地呈现现实中的水体场景。因此还需要对上述模型进行相应的渲染操作，从而得到具有真实感的水体，本发明中采用屏幕空间渲染技术实现水体渲染，从而展示逼真的流体效果。