SPH（光滑粒子流体动力学）流体模拟实现五：PCISPH

SPH（光滑粒子流体动力学）流体模拟实现五：PCISPH

 

我们知道真实的液体是不可压缩的，但我们在计算机中离散的计算流体运动，在一定的时间步长内，用标准的SPH方法求解，在在粒子聚集处容易发生挤压，造成压缩。有两种常用的方法模拟不可压缩性：1.在WCSPH(弱可压缩SPH)中，利用刚性状态方程（EOS）建模压力。2.通过求解压力泊松方程实现不可压缩性。但这两种方法都有很昂贵的计算费用。

文章“Predictive-Corrective Incompressible SPH”中，提出了一种预测矫正的方法，来使粒子达到不可压缩性，其性能上相比传统两种方法，更加的高效。

 

PCISPH模型

SPH概述

在拉格朗日粒子描述下, 控制流体运动的偏微分方程 Navier-Stokes 方程可表示为：

SPH方法的核心思想是以离散化粒子的形式来表征连续的场, 并对场量使用积分近似的方式进行计算。位置在xi粒子i的场量：

第i个粒子的密度计算公式为：

压力场直接从Navier-Stokes方程式推导而得：

PCISPH算法

在PCISPH方法中，速度和位置会及时更新，并估计新的粒子密度。然后，对于每个粒子，计算出参考密度的预测变化，并用于更新压力值，压力值又进入压力的重新计算。此过程一直迭代到收敛，即直到所有粒子密度波动均小于用户定义的阈值η(例如1％)。在完成矫正后，我们再更新速度和位置。详细算法流程图如下：

该算法总体思路如下：

1.计算每个粒子的邻居信息，并记录在邻接表内。

2.计算出了压力之外的所有其他力(黏力，重力)。

3.执行矫正循环：执行1).，2).，3).。

   1).预测所有粒子新的的速度和位置。

   2).预测所有粒子新的密度，以及计算新密度和旧密度之间的差值。

   3).预测所有粒子的压力。

4.完成矫正后，利用矫正的压力计算粒子速度和位置。

 

压强导数

我们的算法是根据预测的密度计算新的压力值，因此我们需要找到它们之间变化的关系。我们的目的是找到一个压强p，该压强以这样一种方式改变粒子的位置，即预测的密度与参考密度相对应。 

对于给定的内核平滑长度h，使用SPH密度求和公式计算时间点t + 1处的密度：

其中，。假定非常小，我们可以用一阶泰勒展开近似：

带入得：

我们令，因此密度增量公式为：

我们假设邻居具有相等的压强且密度为静止密度(根据不可压缩性条件)，则结果是：

PCISPH算法在每次迭代校正时只矫正d流体粒子的压力，通过迭代计算出的压力让粒子之间不至于靠的过近(粒子之间距离太近可以理解为流体可压缩)。在只考虑压力的情况下，根据时间积分可以计算出粒子在该压力作用下的位移：

粒子i在获得粒子j的压力同时，也会对领居粒子j施加反作用力，因此：