【CG物理模拟系列】流体模拟--粒子法之MPS法(理论)

MPS法 

前面的文章里我们讲过SPH曾是为了解决压缩性流体问题而提出的方法，与之相对，这一篇来说说用粒子法处理非压缩性流体的研究方法--Moving Particle Semi-implicit (MPS)法[1]。 MPS法中导入了粒子数密度这一概念，通过隐式方法求解泊松方程，来使得这个密度保持一定。经此处理，流体在整体上便具备了非压缩性。这个方法由于只有一部分使用了隐式方法求解，所以属于Semi-implicit。

非压缩性流体的支配方程

与SPH一样，也用非压缩流体的支配方程(Navier-Stokes方程)来描述流体物质，方程如下.

其中，是流体的速度，是液体密度，为粘度，指外力。

粒子间相互作用模型

MPS法中，梯度，扩散，laplacian等微分算子，是以粒子间相互作用模型为基础进行离散化计算的。 首先加权函数的计算如下.

这里，是粒子间的距离，指代有效半径。这个公式表示，相互作用仅发生在粒子间距小于的情况下。为了保持流体的非压缩性，使用加权函数，按照如下的方式定义粒子数密度。

这里，指代粒子和其临近粒子的位置矢量。另外，我们称粒子数密度的初始值为初期粒子数密度。

由于非压缩性流体的密度是保持一定的。 在MPS法中，通过使各粒子的粒子数密度变成来实现这一特性。

根据粒子间相互作用模型，梯度，扩散，発散，laplacian的离散化如下.

• 梯度：根据粒子的位置得到梯度矢量的演算。
  
  这里，是次元数，是初期粒子数密度.
• 扩散： 
  
• laplacian：把粒子的变量(这里指速度，圧力等物理量)， 根据加权函数的分布，分配到临近的粒子的演算。
  
  其中，的计算如下
  

MPS法中粒子位置的更新

Navier-Stokes方程(质量守恒和动量守恒)的离散化表示如下。

脚标指相应物理量的级数(number of steps)。MPS法中，使用值，計算隐式圧力梯度项。 粘性项和重力项则用显式方法计算。

首先，我们计算粘性项和重力项，获得粒子临时速度，临时位置。

然后通过如下方式离散化粘性项。

这里把粒子位置的粒子数密度表示成，随后，通过解下面的圧力的泊松方程(Poisson)，求得。

和上述粘性项类似，离散化等式右边，可以得到关于的线性系统。

A是对称稀疏矩阵，于是可通过不完全乔勒斯基分解共轭梯度(Incomplete Cholesky Decomposition Conjugate Gradient method)等方法解得。再从得到的中，求解和。

計算算法

MPS法按照如下顺序进行水流模拟。

1. 设定水流的初期条件
2. 显式方法计算临时速度的位置
3. 计算粒子数密度
4. 求解圧力的泊松方程
5. 根据圧力梯度项修正速度的位置
6. 返回步骤2

[1] S. Koshizuka, H. Tamako and Y. Oka, A particle method for incompressible viscous flow with fluid fragmentation, Comput. Fluid Dynamics J., 29(4), 1996.