SPH法理论(1)

SPH法理论

前言

SPH(Smoothed Particle Hydrodynamics)法，中文称其为光滑粒子动力学。最早是为了模拟银河系中天体之间的碰撞与形成等等宇宙物理学的现象而提出的算法[Lucy1977]¹。虽然如此，近些年来SPH法在流体模拟热扩散等等的领域内都被广泛应用了起来。而CG(Computer Graphics)相关的研究者[Desbrun1996]²最早在1996年将SPH法引入计算机模拟“软体”的物理效果。随后[Muller2003]³将SPH法扩展到适用于三维粘性流体的模拟。而本文也将主要基于Muller³在2003年的方法进行介绍与实践。

OpenGL版SSF SPH

Note

本文将涉及到大量数学以及物理知识，限于个人能力如有错误和疑点，请在评论指出。

从粒子角度看世界

既然是光滑粒子动力学，那么自然而然需要从粒子的角度来入手了解这个方法。也就是说，世间万物我们都可以将其看成是由一个一个微小的粒子(质点)组成的。那么对于流体来说，我们并不是将其看作一个整体，比如“水”来看待，而是应该将“水”看作是由无数个微小的粒子而组成的。在这种观点下，我们就需要思考如何计算每一个微小的水粒子的受力以及了解它是如何运动起来的。

光滑粒子动力学

基于上述介绍，我们需要将这种理论先数学模型化，这样才能便于我们的后续分析和研究：
$$f(x_i)=\int f(x_j)\delta (x_i-x_j)d{x}_j(1)$$
此处我们要先思考一个问题，如何计算每一个微小粒子的受力。如果说水是由一个个独立的微小粒子构成的话，那么为了研究粒子$x_i$的受力情况。我们先可以从非常特殊的情况来进行考虑。假设我们先需要获取$x_i$的某一种属性或所受到某种力$f$,即$f(x_i)$,单独来考虑粒子$x_i$是难以计算它的密度或者受到的压强等等的力，但是如果说考虑在粒子$x_i$的周围存在着一群粒子的话，那么$x_i$所受到的力自然而然可以看为是受到了周围粒子$x_j$的作用而产生的。

回到这个数学模型上，$\delta (\cdot )$被称为狄拉克函数，主要有以下的定义：
$$\delta (x)=\{\begin{array}{ccc}\infty & & x=0\\ 0& & x=/0\end{array}$$

$${\int }_{-\infty }^{+\infty }\delta (x)dx=1$$

从Wiki⁴和"小时"⁵所介绍的内容看来，当我们需要独立考察一个质点的密度得到的是无穷大时， 然而其密度在空间中积分却又能得到有限的质量。我们将这样的密度分布可以使用狄拉克函数来表达。再次回到公式$(1)$，我们可以将$\delta (x_i-x_j)$看作是粒子$x_i$和邻近粒子$x_j$之间的关系，而$f(x_j)$则可以表示粒子群$x_j$的某一种属性或者某一种力，那么$f(x_j)\delta (x_i-x_j)$就可以表示为粒子$x_j$对粒子$x_i$产生的某种作用，对其求积分的含义就是粒子$x_i$的某受力$f(x_i)$是由邻近粒子$x_j$所产生某种力$f(x_j)$加和而成，即$\int f(x_j)\delta (x_i-x_j)d{x}_j$。

在分析完理论模型后，我们需要将其离散化这样才可以使得计算机可以去具体进行计算。首先，我们需要将广义的狄拉克函数替换成一个光滑核函数$\omega$。具体来说我们希望可以在考察粒子$x_i$时，能够有一个具体范围内的粒子$x_j$来让我们进行计算。

由此，我们可以将式$(1)$重写为：
$$f(x_i)\approx \int f(x_j)\omega (\frac{\Vert {\mathbf{x}}_i-{\mathbf{x}}_j\Vert }{h})d{p}_j(2)$$

其中$h$是光滑核的半径长度，光滑核函数$\omega$具有以下的性质⁶：

$$\omega (\frac{\Vert {\mathbf{x}}_i-{\mathbf{x}}_j\Vert }{h})=\{\begin{array}{ccc}0& & \frac{\Vert {\mathbf{x}}_i-{\mathbf{x}}_j\Vert }{h}>1\\ \delta ({\mathbf{x}}_i-{\mathbf{x}}_j)& & h\to 0\end{array}$$

$${\int }_{-\infty }^{+\infty }\omega (\frac{\Vert {\mathbf{x}}_i-{\mathbf{x}}_j\Vert }{h})d{x}_j=1$$

当粒子$x_i$和粒子群中某一粒子的距离大于光滑核半径的话，函数值为0。就是为了计算粒子$x_i$的某一属性或者受力情况时我们只考虑来自光滑核内部的粒子群的作用。

在三维空间中，我们的积分变量是体积，因而可以将公式$2$离散成以下形式：
$$\begin{array}{rl}f(x_i)& \approx \int f(x_j)\omega (\frac{\Vert {\mathbf{x}}_i-{\mathbf{x}}_j\Vert }{h})d{x}_j\\ & \approx \sum _{i=1}^N{f}_j\omega (\frac{\Vert {\mathbf{x}}_i-{\mathbf{x}}_j\Vert }{h})V_j\\ & \approx \sum _{i=1}^N\frac{m_j}{{\rho }_j}{f}_j\omega (\frac{\Vert {\mathbf{x}}_i-{\mathbf{x}}_j\Vert }{h})(3)\end{array}$$
式$3$中$V_j$是体积微元, 因为存在$\rho =\frac{m}{V}$，因而可以变换出最终形式。

尽管我们现在可以得到粒子$x_i$的某一属性$f(x_i)$的离散化结果，为了便于后续的计算，此处我们还需要补充如何计算$f(x_i)$的梯度和拉普拉斯算子。

梯度(Gradient)

梯度符号为$\nabla$, 梯度为矢量的空间偏导数，结果依然为矢量。
三维梯度有如下形式:
$$\nabla f(x,y,z)=\left(\frac{\partial f}{\partial x},\frac{\partial f}{\partial y},\frac{\partial f}{\partial z}\right)$$

散度(Divergence)

散度算子仅仅应用于向量场，它衡量在某一点出相应的矢量聚集或者发散程度，测量方向为径向，结果为标量。3维形式的散度算子如下所示：
$$\nabla \cdot \mathbf{u}=\nabla \cdot (u,v,w)=\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}$$
输入是矢量，而输出为标量。类比梯度，散度符号$\nabla \cdot \mathbf{u}$可以理解为梯度$\nabla$与矢量$\vec{u}$的点乘：
$$\nabla \cdot \mathbf{u}=\left(\frac{\partial }{\partial x},\frac{\partial }{\partial y},\frac{\partial }{\partial z}\right)\cdot (u,v,w)=\frac{\partial }{\partial x}u+\frac{\partial }{\partial y}v+\frac{\partial }{\partial z}w$$

拉普拉斯算子(Laplacian)

拉普拉斯算子定义为梯度的散度，符号表示为$\nabla \cdot \nabla$。其中$\nabla \cdot$为散度是标量，后面的$\nabla$为梯度是矢量。拉普拉斯算子即为梯度的散度，为二阶微分算子，三维形式如下：
$$\nabla \cdot \nabla f=\frac{{\partial }^{2}f}{\partial x^2}+\frac{{\partial }^{2}f}{\partial y^2}+\frac{{\partial }^{2}f}{\partial z^2}$$

偏微分方程$\nabla \cdot \nabla f=0$ 被称为拉普拉斯方程。$\nabla \cdot \nabla f=p$被称为泊松方程。

此处补充内容来自于[流体模拟基础]⁷。

物理量$f$的梯度以及拉普拉斯算子

对于公式$3$来说我们将继续做进一步的离散化：
$$\begin{array}{rl}{\nabla }_{i}f(x_i)\approx & {\nabla }_i\int f(x_j)\omega (\frac{\Vert {\mathbf{x}}_i-{\mathbf{x}}_j\Vert }{h})d{x}_j\\ \approx & \int [{\nabla }_{i}f(x_j)]\omega (\frac{\Vert {\mathbf{x}}_i-{\mathbf{x}}_j\Vert }{h})d{x}_j+\int f(x_j)\left[{\nabla }_i\omega (\frac{\Vert {\mathbf{x}}_i-{\mathbf{x}}_j\Vert }{h})\right]d{x}_j\\ \approx & \int f(x_j)\left[{\nabla }_i\omega (\frac{\Vert {\mathbf{x}}_i-{\mathbf{x}}_j\Vert }{h})\right]d{x}_j\\ \approx & \sum _{j=1}^N\frac{m_j}{{\rho }_j}{f}_j\nabla \omega (\frac{\Vert {\mathbf{x}}_i-{\mathbf{x}}_j\Vert }{h})(4)\end{array}$$

求得物理量$f$的梯度后，我们同理也可以得到它的拉普拉斯算子：
$${\nabla }_i^{2}f(x_i)\approx \sum _{j=1}^N\frac{m_j}{{\rho }_j}{f}_j{\nabla }^2\omega (\frac{\Vert {\mathbf{x}}_i-{\mathbf{x}}_j\Vert }{h})(5)$$

但是此处还存在对称性等问题，在后面通过纳韦斯托克斯方程来求粘着项和压力项的时候会简述。

核函数

为了获取粒子群对当前粒子$x_i$影响作用的大小。我们先衡量了粒子群中每一个粒子到粒子$x_i$之间的距离，并将这个距离代入核函数$\omega$中计算出一个数值，这个值在此处可以理解为权重。

对于核函数$\omega$的选取，最常见的是poly6核函数：
$$W_{poly6}(r,d)=\frac{315}{64\pi r^3}{\left(1-\frac{d^2}{r^2}\right)}^{3}0\le d\le r$$
此处的$d$是粒子群中某一粒子$x_j$到粒子$x_i$之间的距离，$r$是核半径。

同时，我们也可以得到他的一阶导数：
$$\nabla W_{poly6}(r,d)=\frac{-945}{32\pi r^5}{\left(1-\frac{d^2}{r^2}\right)}^{2}0\le d\le r$$

以及二阶导数：
$${\nabla }^2{W}_{poly6}(r,d)=\frac{945}{32\pi r^5}\left(1-\frac{d^2}{r^2}\right)\left(\frac{5d^2}{r^2}-1\right)0\le d\le r$$

但是，当我们需要计算核函数的梯度以及拉普拉斯算子时，poly6核函数并不是很好的一个选择。我们为了确保流体的不可压缩性，需要计算流体的压力梯度，而poly6核函数的一阶导数在中心处数值为0这就意味着当两个粒子重合时，由于压力梯度为0而会导致没有力能够使得两个粒子分开，这就违背了流体的不可压缩性。

由于以上原因，普遍在计算梯度或拉普拉斯算子时我们并不会选择poly6光滑核，而是会选择使用spiky核函数：
$$W_{spiky}(r,d)=\frac{15}{\pi r^3}{\left(1-\frac{d}{r}\right)}^{3}0\le d\le r$$
一阶导数：
$$\nabla W_{spiky}(r,d)=\frac{-45}{\pi r^4}{\left(1-\frac{d}{r}\right)}^{2}0\le d\le r$$

二阶导数：
$${\nabla }^2{W}_{spiky}(r,d)=\frac{90}{\pi r^5}\left(1-\frac{d}{r}\right)0\le d\le r$$

同时，我们也提供了光滑核函数的网页绘制版：SPH Kernel 绘制
以及核函数的代码实现(Typescript版)：SPH Kernel 代码

在介绍完核函数之后，我们接下来将会介绍具体的SPH流体模拟算法以及实现。

纳维-斯托克斯方程(Navier-Stokes Equation)

1827年由Navier提出了粘性流体的运动方程，而后在1845年由Stokes独立提出粘性系数为一常数的形式。总的来说，纳维-斯托克斯方程描述了粘性不可压缩流体动量守恒的运动方程，简称N-S方程⁸。而N-S方程又分为两部分，一个是动量方程(式(6))，一个是连续方程(式(14))。先从动量方程的介绍开始：
$$\rho \left(\frac{\partial v}{\partial t}+v\cdot \nabla v\right)=-\nabla p+\nabla \cdot T+f(6)$$

$\rho$为密度，$v$是速度，$t$是时间，$p$是流体压强，$T$为粘滞力，$f$是外力。

动量方程核心其实就是来源于牛顿第二定律：
$$\mathbf{F}=m\mathbf{a}$$
接下来我们开始详细分析纳维-斯托克斯方程的来源。对于一个粒子$p$而言，我们记他的质量为$m$,加速度为$\vec{a}$受力为$\vec{F}$，速度为$\vec{u}\equiv (u(t),v(t),w(t))$，则：
$$\mathbf{a}\equiv \frac{D\mathbf{u}}{Dt}\mathbf{F}=m\frac{D\mathbf{u}}{Dt}$$
此处D是在流体力学中很常用的符号，被称为物质导数(material derivative)⁹。单独定义的原因是u不仅仅是时间的函数，还是空间位置x的函数。由链式法则可以推导出：
$$\frac{D\mathbf{u}}{Dt}=\frac{\partial \mathbf{u}}{\partial t}+\nabla \mathbf{u}\cdot \frac{\partial \mathbf{x}}{\partial t}=\frac{\partial \mathbf{u}}{\partial t}+\nabla \mathbf{u}\cdot \mathbf{u}(7)$$

接下来我们来考虑流体压强，流体压强是流体内部的一种力同时也是维持流体不可压缩的关键⁶。如果说存在一个极小的区域$\Omega$，设其表面积为$A$,压强为$P$ 那么极小的区域$\Omega$所受的压力可以写为:
$${\mathbf{F}}_p=-dP\cdot dA(8)$$
而由牛顿第二定律我们可以将式8写为：
$${\mathbf{F}}_p=m\cdot \mathbf{a}=\rho V\cdot \mathbf{a}=\rho dAdz\cdot \mathbf{a}=-dP\cdot dA(9)$$
此处$z$代表压强差最大方向的长度。由式(9)就可以推出:
$$\mathbf{a}=-\frac{dP}{\rho dz}(10)$$
进而我们就可以得到压力项：
$${\mathbf{F}}_p=m\cdot -\frac{dP}{\rho dz}=-V\cdot \nabla p(11)$$
在推导出压力项后，我们还需要流体之间的粘滞力。对于粘滞力我们可以将其看为当前粒子与周围粒子群之间移动速度的关系。再结合上面补充的拉普拉斯算子可看为梯度的散度也就是可以用来衡量一个场中每点变化量的聚合或者发散程度。那么我们可以用拉普拉斯算子来描述当前粒子与周围粒子之间的移动速度的差距,那么我们可以将粘滞力定义为：
$${\mathbf{F}}_{\nu }=m\cdot \nu {\nabla }^2\mathbf{u}(12)$$
此处$\nu$为粘度系数。

将压力项，粘度项和重力再结合式(7),我们可以得到：
$$m\frac{D\mathbf{u}}{Dt}=m\mathbf{g}-V\cdot \nabla p+m\cdot \nu {\nabla }^2\mathbf{u}$$

$$\frac{D\mathbf{u}}{Dt}=\vec{g}-\frac{1}{\rho }\cdot \nabla p+\nu {\nabla }^2\mathbf{u}$$

$$\frac{\partial \mathbf{u}}{\partial t}+\nabla \mathbf{u}\cdot \mathbf{u}=\mathbf{g}-\frac{1}{\rho }\cdot \nabla p+\nu {\nabla }^2\mathbf{u}(13)$$

由动量守恒的角度可以推导出N-S方程的动量方程形式。但由于还存在物质守恒，我们还需要得到N-S方程的连续方程。具体推导内容可以参考“ 电影工业中的流体模拟（三）”¹⁰：
$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathbf{u})=0(14)$$
对于连续方程来说，首先我们模拟的流体是不可压缩的，因而密度$\rho$是不会随着时间而变化的，即$\frac{\partial \rho }{\partial t}=0$，则可以得到$\nabla \cdot \mathbf{u}=0$。将这个结果代入到式(13)可得：
$$\frac{\partial \mathbf{u}}{\partial t}=\mathbf{g}-\frac{1}{\rho }\cdot \nabla p+\nu {\nabla }^2\mathbf{u}(15)$$
而式(15)正表示了每一个粒子的加速度。到此为止，我们所需要的理论基础已经全部介绍完毕，接下来将结合我们的核函数将加速度中每一项的力进行离散化计算并完成流体的模拟。

流体密度计算

通过式3，我们将物理量$f$设置为密度$\rho$可以得到：
$$\begin{array}{rl}\rho (x_i)& =\sum _{j=1}^N\frac{m_j}{{\rho }_j}{\rho }_j\omega (\frac{\Vert {\mathbf{x}}_i-{\mathbf{x}}_j\Vert }{h})\\ & =\sum _{j=1}^N{m}_j\omega (\frac{\Vert {\mathbf{x}}_i-{\mathbf{x}}_j\Vert }{h})(16)\end{array}$$
在实际尝试中光滑核$\omega$我们选取使用Poly6作为核函数，通过空间探索方法得到$i$点附近核半径为$h$区域内的所有粒子$j$，并根据式16可以计算得出粒子$i$出的密度。

压力项计算

关于压力计算，我们普遍会使用能够体现压力和密度相互关联的状态方程式。理想气体的状态方程式$p=k\rho$，$k$是气体常数。在SPH法相关模拟中应用最多的方法是Desbrun提出的状态方程¹¹：
$$p=k(\rho -{\rho }_0)$$
但是这个方程的缺点是模拟出的流体是可压缩的，如果要计算非压缩性的流体比较常见的是采用泰特方程式¹²。但是同样存在一个问题，如果要应用在实时模拟上那么粒子的位置和速度更新的步长要设定在0.0005左右¹³。这就意味着每模拟1秒我们需要计算2000次左右，离线模拟情况下可以使用。当然还有为了减少步长计算而提出的算法，比如PCISPH¹⁴，此处就不过多展开。因而我们普遍在简单对精度要求不高的情况下都会使用Desbrun提出的方法。

在计算压力项时，如果在计算梯度时使用了式4：
$${\nabla }_{i}f(x_i)\approx \sum _{j=1}^N\frac{m_j}{{\rho }_j}{f}_j\nabla \omega (\frac{\Vert {\mathbf{x}}_i-{\mathbf{x}}_j\Vert }{h})(4)$$
那么就会产生不对称的问题，就是粒子$i$对粒子$j$的力和粒子$j$对粒子$i$的力在计算时不一样。解决方法是对求导函数除以密度函数：
$$\begin{array}{rl}\nabla (\frac{f}{\rho })& =\frac{1}{\rho }\nabla f-\frac{f}{{\rho }^2}\nabla \rho \\ \nabla f& =\rho (\nabla \frac{f}{\rho }+\frac{f}{{\rho }^2}\nabla \rho )(17)\end{array}$$
再将式17离散化后可得到：
$$\begin{array}{rl}\nabla f(x_i)& \approx {\rho }_i\sum _{j=1}^N(\frac{m_j}{{\rho }_j}\frac{f_j}{{\rho }_j}+\frac{f_i}{{\rho }_i^2}\frac{m_i}{{\rho }_j}{\rho }_j)\nabla \omega (\frac{\Vert {\mathbf{x}}_i-{\mathbf{x}}_j\Vert }{h})\\ & \approx {\rho }_i\sum _{j=1}^{N}m(\frac{f_j}{{\rho }_j^2}+\frac{f_i}{{\rho }_i^2})\nabla \omega (\frac{\Vert {\mathbf{x}}_i-{\mathbf{x}}_j\Vert }{h})(18)\end{array}$$
再参考纳维斯托克斯方程（式15）中的压力项：
$$f^{pressure}=-\frac{m}{\rho }\cdot \nabla p(19)$$
我们就可以将式18代入式19后可得：
$$\begin{array}{rl}{f}_i^{pressure}& =-\frac{m_i}{{\rho }_i}{\rho }_i\sum _{j=1}^{N}m(\frac{f_i}{{\rho }_i^2}+\frac{f_j}{{\rho }_j^2})\nabla \omega (\frac{\Vert {\mathbf{x}}_i-{\mathbf{x}}_j\Vert }{h})\\ & =-\sum _{j=1}^N{m}^2(\frac{f_j}{{\rho }_j^2}+\frac{f_i}{{\rho }_i^2})\nabla \omega (\frac{\Vert {\mathbf{x}}_i-{\mathbf{x}}_j\Vert }{h})(20)\end{array}$$

粘着项计算

对于粘着项来说，依然先参考纳维斯托克斯方程（式15）：
$$f^{viscosity}=m\nu {\nabla }^2\mathbf{u}(21)$$
此处计算拉普拉斯算子会产生和计算梯度一样的不对称问题，因而我们可以先对求梯度的函数乘一个密度函数即：
$$\begin{array}{rl}\nabla (\rho f)& =f\nabla \rho +\rho \nabla f\\ \nabla f& =\frac{1}{\rho }\nabla (\rho f)-\frac{f}{\rho }\nabla \rho \\ {\nabla }^{2}f& =\frac{1}{\rho }{\nabla }^2(\rho f)-\frac{f}{\rho }{\nabla }^2\rho (22)\end{array}$$
再将式22离散化后可得到：
$$\begin{array}{rl}{\nabla }^{2}f(x_i)& \approx \frac{1}{{\rho }_j}\sum _{j=1}^N\frac{m_j}{{\rho }_j}{f}_j{\rho }_j{\nabla }^2\omega (\frac{\Vert {\mathbf{x}}_i-{\mathbf{x}}_j\Vert }{h})-\frac{f_i}{{\rho }_j}\sum _{j=1}^N\frac{m_j}{{\rho }_j}{\rho }_j{\nabla }^2\omega (\frac{\Vert {\mathbf{x}}_i-{\mathbf{x}}_j\Vert }{h})\\ & \approx \frac{1}{{\rho }_j}\sum _{j=1}^{N}m(f_j-f_i){\nabla }^2\omega (\frac{\Vert {\mathbf{x}}_i-{\mathbf{x}}_j\Vert }{h})(23)\end{array}$$
继而将式21代入式23中：
$$f_i^{viscosity}=\nu \sum _{j=1}^N{m}^2(\frac{{\mathbf{u}}_j-{\mathbf{u}}_i}{{\rho }_j}){\nabla }^2\omega (\frac{\Vert {\mathbf{x}}_i-{\mathbf{x}}_j\Vert }{h})(24)$$

得出了压力项和粘着项之后，结合纳维斯托克斯方程，我们只需要再将重力以及其他外力(例如：风力)一起计算后便可以得到每个粒子当前时刻所受力的总和。最后再根据当前粒子的位置和速度便可以计算出下一时刻粒子的新的位置和速度，如此反复计算便可以达到模拟流体的效果。

Typescript版简易Demo

本质上来说推导出上面的压力项和粘着项的计算公式后，实际编程实践并不是非常困难。但是关于空间内粒子搜索的算法以及碰撞检测的部分依然需要实现，本文暂不介绍此部分，感兴趣的可以参考一下“FLUID ENGINE DEVELOPMENT”这本书¹⁵。本文实践制作的Demo主要参考了“FLUID ENGINE DEVELOPMENT”¹⁵和”SPH理論”¹³。

在线版本的DEMO：Typescript Version SPH
代码：Typescript Version SPH Code

CUDA版SPH实践

下一篇我会简单介绍如何通过CUDA来实现SPH法和WCSPH法并使用SSF(Screen Space Fluid)方法来实现流体渲染，当然比较常用的方法还有通过MC(Marching Cube)法来实现流体表面重建，可以参考“SPH（光滑粒子流体动力学）流体模拟实现四：各向异性(Anisotropic)表面光滑(1)”¹⁶。

---

1. [Lucy1977]:L. B. Lucy, A Numerical Approach to the Testing of the Fission Hypothesis, The Astronomical Journal, Vol.82, No.12, pp.1013-1024, 1977. ↩︎

2. [Desbrun1996]:M. Desbrun and M.-P. Cani, Smoothed Particles: A New Paradigm for Animating Highly Deformable Bodies, Eurographics Workshop on Computer Animation and Simulation (EGCAS), pp.61-76. 1996. ↩︎

3. [Muller2003]M. Muller, D. Charypar and M. Gross, Particle-based Fluid Simulation for Interactive Applications, In Proc. SCA2003, pp.154-159, 2003. ↩︎ ↩︎

4. Wiki/狄拉克函数 ↩︎

5. 小时物理百科/狄拉克函数 ↩︎

6. 基于SPH的拉格朗日流体模拟 ↩︎ ↩︎

7. 流体模拟Fluid Simulation：流体模拟基础 ↩︎

8. Baidu/纳维-斯托克斯方程 ↩︎

9. 电影工业中的流体模拟（三）- 纳维斯托克斯方程（上） ↩︎

10. 电影工业中的流体模拟（三）- 纳维斯托克斯方程（下） ↩︎

11. [Desbrun1996]:M. Desbrun and M.-P. Cani, Smoothed Particles: A New Paradigm for Animating Highly Deformable Bodies, Eurographics Workshop on Computer Animation and Simulation (EGCAS), pp.61-76. 1996. ↩︎

12. [Becker2007]:M. Becker and M. Teschner, Weakly Compressible SPH for Free Furface Flows, In Proc. SCA2007, pp.209-217, 2007. ↩︎

13. SPH理論 ↩︎ ↩︎

14. [Solenthaler, B]:Solenthaler, B., & Pajarola, R. (2009). Predictive-corrective incompressible SPH. In ACM SIGGRAPH 2009 papers (pp. 1-6). ↩︎

15. [Kim2017]：Fluid Engine Development. Boca Raton: Taylor & Francis, a CRC Press, Taylor & Francis Group，2017. ↩︎ ↩︎

16. SPH（光滑粒子流体动力学）流体模拟实现四：各向异性(Anisotropic)表面光滑(1) ↩︎