SPH算法的理论和实践（1）

前段时间做了一个有关于SPH算法的项目，现在正好抽空把它写出来。SPH（Smoothed Particle Hydrodynamics）是光滑粒子流体动力学方法的意思，说白了就是用粒子模拟流体的流动效果。由于项目当中涉及到利用SPH算法实现多流体混溶模拟，所以下面我会写针对单流体模拟和多流体混溶模拟分别进行介绍，在介绍之前还是先看看做出来的最终效果吧。

效果预览

单流体模拟的效果：

多流体混溶模拟的效果：

针对上述多流体混溶模型，利用marching cube算法进行液面绘制之后的效果：

1单流体理论部分

SPH (Smoothed Particle Hydrodynamics)是光滑粒子流体动力学方法的缩写，是在近20多年来逐步发展起来的一种无网格方法。该方法的基本思想是将连续的流体用相互作用的粒子来描述，各个粒子上承载各种物理量，包括质量、速度等，通过求解粒子的动力学方程和跟踪每个粒子的运动轨道，求得整个系统的力学行为。

1.1粒子受力分析

SPH方法具体的数学公式可以看这篇文章，里面介绍的很详细，我当时也是看这篇文章入门的。不过在这里我还是简单介绍一下里面的几个基本公式：

SPH算法的基本思想，就是将连续的流体想象成一个个相互作用的粒子，这些粒子相互影响，共同形成了复杂的流体运动。

所以，该算法最关键的就是如何求解粒子的运动？我们通过对粒子进行受力分析，利用牛顿第二定理，只要知道了粒子的受力，粒子的加速度自然就知道了，这样就可以跟踪粒子的运动，进而模拟出整个流体的动态效果。那如何求解粒子的力F 呢？

根据流体动力学，作用在粒子上的力由三部分组成：

其中，表示外力，一般就是重力

注意：这里的力F的量纲发生了变化，正常情况下，F=m*g。但在SPH算法里，流体的质量是由流体单元的密度决定的，所以一般用密度代替质量，后面的分析都是用这个量纲的“作用力”

表示压力，是由流体内部的压力差产生的作用力，试想一下在水管中流动的液体，进水口区域的压力一定会比出水口区域大，所以液体才会源源不断的流动，数值上，它等于压力场的梯度，方向由压力高的区域指向压力低的区域

表示粘滞力，是由粒子之间的速度差引起的，设想在流动的液体内部，快速流动的部分会施加类似于剪切力的作用力到速度慢的部分，这个力的大小跟流体的粘度系数μ以及速度差有关

所以，最后得到粒子的受力公式

加速度形式：

至此，我们基本搞清楚了粒子的运动学计算方法。下面介绍SPH算法的关键部分，如何通过光滑核函数求解上述公式

1.2光滑核函数

和其他流体力学中的数学方法类似，SPH算法同样涉及到“光滑核”的概念，可以这样理解这个概念，粒子的属性都会“扩散”到周围，并且随着距离的增加影响逐渐变小，这种随着距离而衰减的函数被称为“光滑核”函数，最大影响半径为“光滑核半径”。

反过来不难理解，尽管我们将流体视为一个个分散的粒子，但流体毕竟是连续充满整个空间的，流体中每个位置参与运算的值都是由周围一组粒子累加起来的。

设想流体中某点 r⃗ r → \vec r（此处不一定有粒子）,在光滑核半径h范围内有数个粒子，位置分别是 r⃗ 0,r⃗ 1,r⃗ 2,...,r⃗ j r → 0 , r → 1 , r → 2 , . . . , r → j \vec r_0,\vec r_1,\vec r_2,...,\vec r_j，则该处某项属性A的累加公式为：

其中 A⃗ j A → j \vec A_j是要累加的某种属性， mj m j m_j和 ρj ρ j \rho_j是周围粒子的质量和密度， r⃗ r → \vec r 是该粒子的位置， h h h是光滑核半径。函数W" role="presentation" style="position: relative;">WWW就是光滑核函数。
光滑核函数两个重要属性，首先一定是偶函数，也就是 W(−r)=W(r) W ( − r ) = W ( r ) W(−r)=W(r)，第二，是“规整函数”，也就是 ∫W(r)dr=1 ∫ W ( r ) d r = 1 ∫W(r)dr=1

SPH推导过程
假设流体中的一个粒子 r⃗ i r → i \vec r_i,假设它的压力为 ρ(ri) ρ ( r i ) \rho(r_i)，密度为 p(ri) p ( r i ) p(r_i),速度为 u⃗ (ri) u → ( r i ) \vec u(r_i),name我们可以根据上一节公式，得到该粒子的加速度 a⃗ (ri) a → ( r i ) \vec a(r_i)：

a⃗ (ri)=g⃗ −∇p(ri)ρ(ri)+μ∇2u(ri)ρ(ri) a → ( r i ) = g → − ∇ p ( r i ) ρ ( r i ) + μ ∇ 2 u ( r i ) ρ ( r i )

\vec a(r_i)=\vec g - \frac{∇p(r_i)}{ρ(r_i)}+\frac{μ∇^2u(r_i)}{ρ(r_i)}
对于SPH算法来说，基本流程就是这样，根据光滑核函数逐个推出流体中某点的密度，压力，速度相关的累加函数，进而推导出此处的加速度，从而模拟流体的运动趋势。

下面我们直接利用光滑核函数求解某点的密度，压力，速度的公式，具体的数学推导可以看这里。

密度
根据光滑核函数，用密度 ρ ρ \rho代替 A A A,可以得到

ρ(ri)=∑jρjmjρjW(r→i−r→j,h)=∑jmjW(r→i−r→j,h)" role="presentation">ρ(ri)=∑jρjmjρjW(r⃗ i−r⃗ j,h)=∑jmjW(r⃗ i−r⃗ j,h)ρ(ri)=∑jρjmjρjW(r→i−r→j,h)=∑jmjW(r→i−r→j,h)

\rho(r_i)=\sum_j\frac{ρ_j}{m_j}ρ_jW(\vec r_i−\vec r_j,h)=\sum_jm_jW(\vec r_i−\vec r_j,h)

利用Poly6函数得到

ρ(ri)=m31564πh9∑j(h2−|r⃗ i−r⃗ j|)3 ρ ( r i ) = m 315 64 π h 9 ∑ j ( h 2 − | r → i − r → j | ) 3

\rho(r_i)=m\frac{315}{64\pi h^9}\sum_j (h^2 - |\vec r_i - \vec r_j|)^3

压强
对于单个粒子的压力 p p p可以利用理想气体状态方程计算

p=K(ρ−ρo)" role="presentation">p=K(ρ−ρo)p=K(ρ−ρo)

p = K(\rho - \rho_o)
其中 ρ0 ρ 0 \rho_0是流体的静态密度， K K K是和流体相关的常数，只跟温度有关。

压力
根据光滑核函数，用压力p" role="presentation" style="position: relative;">ppp代替 A A A,可以得到

Fipressure=−∇p(r→i)=−∑jpjmjρj∇W(r→i−r→j,h)" role="presentation">Fpressurei=−∇p(r⃗ i)=−∑jpjmjρj∇W(r⃗ i−r⃗ j,h)Fipressure=−∇p(r→i)=−∑jpjmjρj∇W(r→i−r→j,h)

F_i^{pressure}=-∇p(\vec r_i)=-\sum_j p_j \frac {m_j}{\rho_j} ∇W(\vec r_i - \vec r_j,h)
利用Spiky函数，得到

apressurei=−∇p(r⃗ i)ρ0=m45πh6∑j(pi+pj2ρIρj(h−r)2r⃗ i−r⃗ jr) a i p r e s s u r e = − ∇ p ( r → i ) ρ 0 = m 45 π h 6 ∑ j ( p i + p j 2 ρ I ρ j ( h − r ) 2 r → i − r → j r )

a_i^{pressure} = -\frac{∇p(\vec r_i)}{\rho_0}=m\frac{45}{\pi h^6}\sum_j (\frac{p_i+p_j}{2\rho_I\rho_j}(h-r)^2\frac{\vec r_i - \vec r_j}{r})
其中 r=|r⃗ i−r⃗ j| r = | r → i − r → j | r = |\vec r_i - \vec r_j|

粘滞力
根据光滑核函数，可以得到

Fviscosityi=μ∇2u⃗ (ri)=μ∑ju⃗ jmjρj∇W(r⃗ i−r⃗ j,h) F i v i s c o s i t y = μ ∇ 2 u → ( r i ) = μ ∑ j u → j m j ρ j ∇ W ( r → i − r → j , h )

F_i^{viscosity}=\mu∇^2\vec u(r_i)=\mu \sum_j \vec u_j \frac {m_j}{\rho_j} ∇W(\vec r_i - \vec r_j,h)
根据viscosity函数，得到

aviscosityi=Fviscosityiρ0=μ∇2u⃗ (ri)ρ0=mμ45πh6∑j(u⃗ j−u⃗ iρIρj(h−r)) a i v i s c o s i t y = F i v i s c o s i t y ρ 0 = μ ∇ 2 u → ( r i ) ρ 0 = m μ 45 π h 6 ∑ j ( u → j − u → i ρ I ρ j ( h − r ) )

a_i^{viscosity} = \frac{F_i^{viscosity}}{\rho_0}=\frac{\mu∇^2\vec u(r_i)}{\rho_0}=m\mu\frac{45}{\pi h^6}\sum_j (\frac{\vec u_j - \vec u_i}{\rho_I\rho_j}(h-r))
其中 r=|r⃗ i−r⃗ j| r = | r → i − r → j | r = |\vec r_i - \vec r_j|

粒子的运动方程
将上述压力和粘滞力带入中，得到

a⃗ (ri)=g⃗ +m45πh6∑j(pi+pj2ρIρj(h−r)2r⃗ i−r⃗ jr)+mμ45πh6∑j(u⃗ j−u⃗ iρIρj(h−r)) a → ( r i ) = g → + m 45 π h 6 ∑ j ( p i + p j 2 ρ I ρ j ( h − r ) 2 r → i − r → j r ) + m μ 45 π h 6 ∑ j ( u → j − u → i ρ I ρ j ( h − r ) )

\vec a(r_i)=\vec g + m\frac{45}{\pi h^6}\sum_j (\frac{p_i+p_j}{2\rho_I\rho_j}(h-r)^2\frac{\vec r_i - \vec r_j}{r}) + m\mu\frac{45}{\pi h^6}\sum_j (\frac{\vec u_j - \vec u_i}{\rho_I\rho_j}(h-r))
其中 r=|r⃗ i−r⃗ j| r = | r → i − r → j | r = |\vec r_i - \vec r_j|

OK，我们终于推导除了粒子加速度 a⃗ (ri) a → ( r i ) \vec a(r_i)的计算公式。那如何让这些公式在代码里跑起来呢？不用担心，这就是下一章我们要解决的问题了。