流体模拟：Smoothed Particle Hydrodynamics
SPH，全称为Smoothed Particle Hydrodynamics，最初提出于天体物理学领域，然后被广泛的应用到计算流体力学领域，成为基于拉格朗日粒子模拟方法的典型代表。实际上，目前除了流体，还有弹性固定、刚性固体等的物理模拟也有不少采用了SPH的方法。SPH是一种基于光滑粒子核的物理模型，它将模拟的对象离散成一个个粒子，然后以光滑核将粒子之间联系起来，显然这是一种基于拉格朗日视角的模拟方法，相对于欧拉视角的模拟方法，它比较简单、速度较快。
1.光滑粒子动力学（SPH）
SPH 通常被理解为空间场量和空间微分算子(如梯度、散度、旋度等)离散化的一种方法。为了理解这个想法我们先介绍狄拉克δ函数：
\[\delta(\mathbf{r})= \begin{cases}\infty & \text { if } \quad \mathbf{r}=\mathbf{0} \\ 0 & \text { otherwise }\end{cases} \tag{1} \label{1}
\]
且同时满足
\[\int \delta(\mathbf{r}) d v=1
\]
狄拉克δ函数可理解为方差\(\sigma^2 \rightarrow 0\)的标准高斯分布。
\[\mathcal{N}\left(\mathbf{x} ; \mu, \sigma^{2}\right), \lim _{\sigma \rightarrow 0} \mathcal{N}\left(\mathbf{x} ; 0, \sigma^{2}\right)=\delta(\mathbf{x})
\]
我们使用狄拉克δ函数作为离散化的基础，我们设空间中处于位置x的粒子的物理量为\(A(x)\)，那么\(A(x)\)在空间可以近似表示为：
\[A(\mathbf{x})=(A * \delta)(\mathbf{x})=\int A\left(\mathbf{x}^{\prime}\right) \delta\left(\mathbf{x}-\mathbf{x}^{\prime}\right) d v^{\prime} \tag{2} \label{2}
\]
\(dv'\)代表的是对变量\(x'\)的体积微元，我们用一个光滑的内核函数（光滑因子）来替换狄拉克δ函数：
\[\begin{aligned}
A(\mathbf{x}) & \approx(A * W)(\mathbf{x}) \\
&=\int A\left(\mathbf{x}^{\prime}\right) W\left(\mathbf{x}-\mathbf{x}^{\prime}, h\right) d v^{\prime}
\end{aligned}
\tag{3} \label{3}
\]
其中\(h\)是光滑核的半径长度，光滑核函数\(W\)应该具有以下的性质：
-
归一化条件：\(\int_{\mathbb{R}^{d}} W\left(\mathbf{r}^{\prime}, h\right) d v^{\prime}=1\)
-
狄拉克δ函数条件：当核半径\(h\)趋于零时，核函数收敛到狄拉克δ函数
\[\lim _{h^{\prime} \rightarrow 0} W\left(\mathbf{r}, h^{\prime}\right)=\delta(\mathbf{r})
\]
-
正性条件：\(W(r,h) \ge 0\),即超出核半径之外的定义域，函数值均取零
-
对称条件：\(W(-r,h)=W(r,h)\)，保证了连续近似的一阶一致性。
-
函数\(W\)是光滑连续的，因而是可微的
-
函数\(W\)应该是二阶连续可微，以保证二阶微分方程的一致离散化
通过公式\((3)\)，我们可以根据周围邻域的函数值近似获取给定的点的函数值，这些函数值可以是密度、压强等流体的属性值。公式\((3)\)给出的是一个积分形式，在二维空间中的积分变量是面积，在三维空间中的积分变量是体积，我们采用黎曼和来离散地计算公式\((3)\)：
\[\begin{aligned}
(A * W)\left(\mathbf{x}_{i}\right) &=\int \frac{A\left(\mathbf{x}^{\prime}\right)}{\rho\left(\mathbf{x}^{\prime}\right)} W\left(\mathbf{x}-\mathbf{x}^{\prime}, h\right) \underbrace{\rho\left(\mathbf{x}^{\prime}\right) d v^{\prime}}_{d m^{\prime}} \\
& \approx \sum_{j \in \mathcal{F}} A_{j} \frac{m_{j}}{\rho_{j}} W\left(\mathbf{x}_{i}-\mathbf{x}_{j}, h\right)=:\left\langle A\left(\mathbf{x}_{i}\right)\right\rangle
\end{aligned}
\tag{4} \label{4}
\]
其中\(dv'\)是我们取的体积微元，\(\mathcal{F}\)表示所有点集，,\(W_{i j}=W(x_i-x_j,h)\),又因\(\rho=\frac{m}{v}\),所以得到SPH方法的基本公式\((4)\),该公式可以用于近似问题域中任何位置处的任何连续的场量。同时因为：
\[\begin{aligned}
\left\langle A\left(\mathbf{x}_{i}\right)\right\rangle=& A_{i} \sum_{j} \frac{m_{j}}{\rho_{j}} W_{i j}+\\
&\left.\nabla A\right|_{\mathbf{x}_{i}} \cdot \sum_{j} \frac{m_{j}}{\rho_{j}}\left(\mathbf{x}_{j}-\mathbf{x}_{i}\right) W_{i j}+\mathcal{O}\left(\left(\mathbf{x}_{j}-\mathbf{x}_{i}\right)^{2}\right)
\end{aligned}
\]
如果下式成立，
\[\sum_{j} \frac{m_{j}}{\rho_{j}} W_{i j}=1 \quad \text { and } \quad \sum_{j} \frac{m_{j}}{\rho_{j}}\left(\mathbf{x}_{j}-\mathbf{x}_{i}\right) W_{i j}=\mathbf{0}
\]
那么离散化的结果只有一阶精度。
1.1质量密度估计
粒子并不需要携带质量密度信息，我们使用SPH基本方程\((4)\)根据位置来估计质量密度
\[\rho_{i}=\sum_{j} m_{j} W_{i j} \tag{5} \label{5}
\]
1.2 微分算子离散化
为了物理公式进行数值模拟，我们还需要对空间微分算子进行离散化操作。根据SPH基本方程\((4)\)，由于\(A_j \frac{m_j}{\rho_j}\),不依赖于积分变量，可以视为常量积分，因此微分算子只影响SPH核函数，因此我们可以得到SPH公式的梯度:
\[\nabla A_{i} \approx \sum_{j} A_{j} \frac{m_{j}}{\rho_{j}} \nabla W_{i j}
\]
甚至一些更为复杂的空间微分算子：
\[\begin{aligned}
\nabla \mathbf{A}_{i} & \approx \sum_{j} \frac{m_{j}}{\rho_{j}} \mathbf{A}_{j} \otimes \nabla W_{i j} \\
\nabla \cdot \mathbf{A}_{i} & \approx \sum_{j} \frac{m_{j}}{\rho_{j}} \mathbf{A}_{j} \cdot \nabla W_{i j} \\
\nabla \times \mathbf{A}_{i} & \approx-\sum_{j} \frac{m_{j}}{\rho_{j}} \mathbf{A}_{j} \times \nabla W_{i j}
\end{aligned}
\]
\(a\otimes b=ab^T\)代表并矢积（二元乘积），但是上述方程“直接”导数不好的近似结果和不稳定的模拟。所以我们提出了以下使用最广泛的SPH一阶微分算子离散公式。
差分方程：
根据泰勒级数对梯度误差进行分析，如下满足以下第一个条件（都满足），梯度近似值只有0阶（1阶）精确：
\[\langle\nabla 1\rangle=\sum_{j} \frac{m_{j}}{\rho_{j}} \nabla_{i} W_{i j}=\mathbf{0} \quad \text { and } \quad \sum_{j} \frac{m_{j}}{\rho_{j}}\left(\mathbf{x}_{j}-\mathbf{x}_{i}\right) \otimes \nabla_{i} W_{i j}=1
\]
为了恢复0阶精确，我们需要简单地减去泰勒级数的一个误差项，从而得到一个梯度的差分公式：
\[\nabla A_{i} \approx\left\langle\nabla A_{i}\right\rangle-A_{i}\langle\nabla 1\rangle=\sum_{i} \frac{m_{j}}{\rho_{j}}\left(A_{j}-A_{i}\right) \nabla_{i} W_{i j}
\tag{6} \label{6}
\]
对称方程：
根据流体力学系统，用拉格朗日离散方法和密度估计，我们可以推断出一个压力（梯度）的方程：
\[\begin{aligned}
\nabla A_{i} & \approx \rho_{i}\left(\frac{A_{i}}{\rho_{i}^{2}}\langle\nabla \rho\rangle+\left\langle\nabla\left(\frac{A_{i}}{\rho_{i}}\right)\right\rangle\right) \\
&=\rho_{i} \sum_{j} m_{j}\left(\frac{A_{i}}{\rho_{i}^{2}}+\frac{A_{j}}{\rho_{j}^{2}}\right) \nabla_{i} W_{i j}
\end{aligned}
\tag{7} \label{7}
\]
可以看出，并不满足之前的两个约束条件，因此，它不能够准确地重现常数线性梯度方程；但是，使用这种特殊的梯度估计的离散物理力是保守的线性和角动量，这是鲁棒模拟的一个基本准则。
用泰勒级数展开式推导判据，对称梯度方程的常数误差取决于0偏离多少：
\[\sum_{j} m_{j}\left(\frac{1}{\rho_{i}^{2}}+\frac{1}{\rho_{j}^{2}}\right) \nabla_{i} W_{i j} \approx \mathbf{0}
\]
对称方程“关心”的是粒子的有序性和试图重新派刘的粒子结构的物理力知道上式满足。这与使用差分方程表示的相反。
因此，差分方程的更为精确的梯度估计是用动量的损失来实现的，因此会导致不稳定的模拟。
在使用例如物理力等直接影响粒子轨迹的量，我们建议使用对称方程进行梯度估计。
1.3 拉普拉斯算子离散化
与直接一阶求导相等，拉普拉斯算子可以直接用下述公式描述：
\[\nabla^{2} A_{i} \approx \sum_{j} \frac{m_{j}}{\rho_{j}} A_{j} \nabla_{i}^{2} W_{i j}
\]
然而，这又导致了对二阶的以一个糟糕的估计，因此Brookshaw提出了一种改进的拉普拉斯算子：
\[\nabla^{2} A_{i} \approx-\sum_{j} \frac{m_{j}}{\rho_{j}} A_{i j} \frac{2\left\|\nabla_{i} W_{i j}\right\|}{\left\|\mathbf{r}_{i j}\right\|}
\]
上式公式有效仅使用核函数得一阶导数，并使用类似有限差分的差分运算（即除以粒子距离）来实现二阶导数。
用类比的方法是实现了矢量场量的二阶导数：
\[\begin{aligned}
\nabla^{2} \mathbf{A}_{i} &=-\sum_{j} \frac{m_{j}}{\rho_{j}} \mathbf{A}_{i j} \frac{2\left\|\nabla_{i} W_{i j}\right\|}{\left\|\mathbf{r}_{i j}\right\|} \\
\nabla\left(\nabla \cdot \mathbf{A}_{i}\right) &=\sum_{j} \frac{m_{j}}{\rho_{j}}\left[(d+2)\left(\mathbf{A}_{i j} \cdot \tilde{\mathbf{r}}_{i j}\right) \tilde{\mathbf{r}}_{i j}-\mathbf{A}_{i j}\right] \frac{\left\|\nabla_{i} W_{i j}\right\|}{\left\|\mathbf{r}_{i j}\right\|}
\end{aligned}
\]
\(d\)表示空间维度，\(\mathbf{r}\)表示粒子距离的标准化向量。
但是上式第一个公式对一些力（如粘性力）并不是动量守恒的；不过我们通过上述两个公式相加得到
\[\sum_{j} \frac{m_{j}}{\rho_{j}}\left(\mathbf{A}_{i j} \cdot \tilde{\mathbf{r}}_{i j}\right) \tilde{\mathbf{r}}_{i j} \frac{\left\|\nabla_{i} W_{i j}\right\|}{\left\|\mathbf{r}_{i j}\right\|} \approx \frac{\nabla\left(\nabla \cdot \mathbf{A}_{i}\right)}{d+2}-\frac{\nabla^{2} \mathbf{A}_{i}}{2(d+2)}
\]
这个恒等于一个作用是：在无散度的向量场\(\nabla \cdot A=0\)，拉普拉斯算子可以使用，如下式：
\[\nabla^{2} \mathbf{A}_{i} \approx 2(d+2) \sum_{j} \frac{m_{j}}{\rho_{j}} \frac{\mathbf{A}_{i j} \cdot \mathbf{r}_{i j}}{\left\|\mathbf{r}_{i j}\right\|^{2}} \nabla_{i} W_{i j} \tag{8} \label{8}
\]
这样可以使物理力导出的力符合动量守恒。因此，我们使用上式进行对拉普拉斯算子的离散化。
1.4 核函数
1.4.1 三次样条插值函数
光滑核的一个典型选择是三次样条插值函数：
\[W(\mathbf{r}, h)=\sigma_{d} \begin{cases}6\left(q^{3}-q^{2}\right)+1 & \text { for } 0 \leq q \leq \frac{1}{2} \\ 2(1-q)^{3} & \text { for } \frac{1}{2}<q \leq 1 \\ 0 & \text { otherwise }\end{cases}
\]
其中\(q=\frac{1}{h} ||r||\),各维度的正规化因子为\(\sigma_1=\frac{4}{3h},\sigma_2=\frac{40}{7\pi h^2},\sigma_3=\frac{8}{\pi h^3}\)
下面介绍下常用于流体模拟的核函数：
1.4.2 Poly6核函数
\[W_{p o l y 6}(r, h)=\frac{315}{64 \pi h^{3}} \begin{cases}\left(1-\frac{r^{2}}{h^{2}}\right)^{3} & 0 \leq r \leq h \\ 0 & \text { otherwise }\end{cases}
\]
其梯度计算公式为：
\[\nabla W_{p o l y 6}(r, h)=\frac{-945 r}{32 \pi h^{5}}\left(1-\frac{r^{2}}{h^{2}}\right)^{2}, 0 \leq r \leq h
\]
二阶导数（拉普拉斯算子）为：
\[\nabla^{2} W_{p o l y 6}(r, h)=\frac{945}{32 \pi h^{5}}\left(1-\frac{r^{2}}{h^{2}}\right)\left(\frac{5 r^{2}}{h^{2}}-1\right)
\]
当我们需要求取核函数的梯度向量核拉普拉斯算子的时候，我们并不会采用poly6核函数，这是因为虽然其原函数很光滑，但其一阶导数和二阶导数的性质却不是很好，因而不会使用poly6的导数。因此对于梯度向量和拉普拉斯的求取，我们用spiky核函数取而代之。
1.4.3 spkiy核函数
\[W_{spkiy}(r, h)=\frac{15}{\pi h^{3}} \begin{cases}\left(1-\frac{r}{h}\right)^{6} & 0 \leq r \leq h \\ 0 & \text { otherwise }\end{cases}
\]
下面的图是关于poly6核函数与spiky核函数的函数曲线、一阶导数曲线和二阶导数曲线，左图是poly6函数的，右图是spiky函数的，实心曲线为函数自身的曲线，点线是函数的二阶导数曲线，而虚线则是函数的一阶导数曲线。可以看到poly6的一阶导数和二阶导数呈现出一个波动的形态，在中心处一阶导数甚至降为零。我们需要计算流体的压力梯度来确保流体的不可压缩性，如果直接采用poly6的一阶导数来计算流体的压力梯度力，那么当两个粒子重合时，其压力梯度力为零，不存在一个力使它们分开，从而违背了流体的不可压缩性。二阶导数同理，其甚至在中心处取值为负，二阶导数将用于拉普拉斯算子的计算，而流体的粘性力计算将用到拉普拉斯算子
spiky函数的一阶导数为：
\[\nabla W_{\text {spiky }}(r, h)=-\frac{45}{\pi h^{4}}\left(1-\frac{r}{h}\right)^{2}\frac{r}{|r|}
\]
spiky函数的二阶导数为：
\[\nabla^{2} W_{\text {spliky }}(r, h)=-\frac{90}{\pi h^{6}} \frac{1}{r}(h-r)(h-2 r)
\]
因此poly6核函数常用于描述密度\(\rho\),而skipy核函数常用于描述压力\(p\)
1.4.4 viscosity核函数
\[W_{\text {viscosity }}(\mathbf{r}, h)=\frac{15}{2 \pi h^{3}} \begin{cases}-\frac{r^{3}}{2 h^{3}}+\frac{r^{2}}{h^{2}}+\frac{h}{2 r}-1 & 0 \leq r \leq h \\ 0, & \text { otherwise }\end{cases}
\]
viscosity核函数常用于粘滞力\(\mu\)计算，粘滞现象主要有摩擦力引起，因此会减少流体的动能（转化为热能）。因此，粘滞性项的计算应只对流体速度场有平滑作用，如果采用poly6核函数则计算产生的粘滞力并不一定会具备这个性质。
viscosity核函数的一阶导数（梯度）为：
\[\nabla W_{\text {viscosity }}(r, h)=\frac{15}{2 \pi h^{3}} r\left(-\frac{3 r}{2 h^{3}}+\frac{2}{h^{2}}-\frac{h}{2 r^{3}}\right)
\]
viscosity核函数的二阶导数（拉普拉斯算子）为：
\[\nabla^{2} W_{v i s c o s i t y}(r, h)=\frac{45}{r h^{6}}(h-r)
\]
具体如下图：
2.SPH流体模拟算法
2.1 连续性方程(Continuity Equation)
连续性定理：物质是由离散的分子、原子构成，但在宏观上可以被理想化为一个连续体，它的任一部分总能被划分为更小的部分。连续介质理论忽略微观的影响，着眼于宏观物理现象。
\[\frac{\partial \rho}{\partial t}+\nabla \cdot(\rho \vec{v})=0
\tag{9} \label{9}
\]
方程\((9)\)描述了流体的密度随着时间变化速率，我们模拟的是不可压缩的流体，因而密度守恒，所以\(\frac{\partial \rho}{\partial t}=0\),即\(\nabla \cdot \vec{v}=0\)(不可压缩约束条件)，流体的速度场无散度。
\(\frac{D}{Dt}\)为物质导数(Material Derivatvie),物质导数描述了一个场量在某个质点处关于时间的导数(随时间变化率)。物质导数的显式形式取决于问题域空间类型，在欧拉流体问题域物质导数描述空间某固定点处的场量关于时间的变化率：
\[\frac{D A^{E}}{D t}=\frac{\partial A^{E}}{\partial t}+v \cdot \nabla_{x} A^{E}
\]
其中\(v \cdot \nabla_{x} A^{E}\)就是欧拉流体中的平流项（Convection Term）。
在拉格朗日流体问题域中物质导数描述了空间中某粒子的场量关于时间的变化率：
\[\frac{D A^{L}}{D t}=\frac{\partial A^{L}}{\partial t}
\]
显然拉格朗日流体中不存在平流项。
2.2 线动量守恒定律(Consevation Law of Linear Momentum)
表明粒子动量的变化率等于作用于粒子的所有内外力之和：
\[\rho \frac{D^{2} x}{D t^{2}}=\nabla \cdot T+f_{e x t}
\]
其中\(f_{ext}\)为体积力（Body Force）,在此处指的是单位体积力而不是物质整体的体积力。\(T\)为应力张力（Steess Tensor），由物质本身的特性决定，它与物质自身的本构关系有关，不可压缩流体的经典本构关系为\(T=-pE+\mu (\nabla v+\nabla v^T)\),根据这个关系可以得到不可压缩流体的Navier-Stokes方程\((10)\)
\[\rho \frac{D v}{D t}=-\nabla p+\mu \nabla^{2} v+f_{\text {ext }}
\tag{10} \label{10}
\]
N-S方程是流体力学中关于不可压缩流体的经典方程，表示了流体受到压力、粘性力以及其他外力的作用。如果不可压缩性不严格(微可压)，则压强\(p\)可以通过状态方程(State Eq.)施加，压力与流体体积构成本构关系\(p=p(\rho)\),状态方程有多种，如理想气体状态方程（Ideal Gas State Eq.）的变形：\(p=k(\rho - \rho _0)\)或者\(\rho = k(\frac{\rho}{\rho_0}-1)\)就是把密度偏差和压力联系起来，并通过刚度系数\(k\)来线性放缩得到压强。
2.3 不可压缩流体的Navier-Stokes方程
上式左边\((10)\)的\(\frac {D\vec{v}}{Dt}\)是速度关于时间的物质导数，即\(\frac {D\vec{v}}{Dt}=\frac{\partial \vec{v}}{\partial t}+\vec{v} \cdot \nabla \vec{v}\)，在我们基于粒子的流体模拟算法中，质量守恒满足即\(\vec{v} \cdot \nabla \vec{v}=0\),所以粒子速度场关于时间的物质导数就是\(\frac{\partial \vec{v}}{\partial t}\);现在体积力项为\(f_{ext}=\rho \vec{g}\)(只考虑重力)，\(\mu\)是流体的粘度系数，同时定义运动粘度为\(v=\frac{\mu}{\rho}\)，\(\rho\)是流体密度，\(p\)是流体压强，可得到下述公式\((11)\)：
\[\frac{D \vec{v}}{D t}=-\frac{1}{\rho} \nabla p+v \nabla^{2} \vec{v}+\vec{g}
\tag{11} \label{11}
\]
每一个粒子的加速度为：\(a=-\frac{1}{\rho} \nabla p+v \nabla^{2} \vec{v}+\vec{g}\);方程右边分别是流体的压力项，粘性力项，以及体积力项。
所以最终的Navier-Stokes方程即为：
\[\begin{eqnarray*}
\frac{D \vec{v}}{D t}=-\frac{1}{\rho} \nabla p+\mu \nabla^{2} \vec{v}+\vec{g} \tag{10}\\
\frac{\partial \rho}{\partial t}+\nabla \cdot(\rho \vec{v})=0 \tag{11}
\end{eqnarray*}
\]
2.3 算法过程
SPH模拟流体实际上是通过SPH方法求NS方程的数值解，为了确保解的唯一性，需要有初值条件，它可能包括流体的初始形状(每个粒子的初始位置)、每个粒子的初始速度、位置的约束条件(空间边界)或速度的约束条件等。因此SPH模拟流体即构成一个初值-边界值混合问题(Mixed Initial-Boundary Value Problem)。
在模拟的过程中，我们需要分别计算这几项。一个完整的SPH流体求解器计算流程如下所示:
1.Measure the density with particles’ current locations
2.Compute the gravity and other extra forces
3.Compute the viscosity force
4.Compute the pressure based on the density
5.Compute the pressure gradient force
6.Perform time integration
接下来按照这个步骤一一展开：
1丶计算粒子的密度
首先我们要计算流体粒子的密度值，因为后面粘性力和压强梯度力的计算将会用到粒子的密度。粒子的额密度同样是采用光滑核对周围粒子进行加权计算，将密度函数代入公式\((4)\):
\[\begin{aligned}
f(x) & \approx \Sigma_{i=1}^{N} \frac{m_{i}}{\rho_{i}} \rho_{i} w\left(\left\|x-x_{i}\right\| / h_{i}\right) \\
& \approx \Sigma_{i=1}^{N} m_{i} w\left(\left\|x-x_{i}\right\| / h_{i}\right)
\end{aligned}
\]
原公式中的密度项被消去了，所以我们直接根据粒子的位置和质量计算每个粒子的密度，在这里我们的每个粒子质量都相同。
2丶计算外部作用力
这里指的外部的作用力通常是体积力，即间接地作用在流体上而非通过直接接触产生的作用力。常见的体积力有重力和风力，这些体积力我们直接根据需要指定其加速度的值，然后采用牛顿定律叠加到粒子的力场上
3丶计算流体粘性力
流体的粘性力是流体内部的一种阻力，粘性力的计算需要计算流体速度场的拉普拉斯算子，这是因为拉普拉斯算子衡量了给定位置的物理量与周围领域的物理量的差距值：
\[f_{v}=m \mu \nabla^{2} \vec{v}
\]
在SPH的算法中，其同样是通过光滑核函数计算粘性力，注意为了避免常量函数的拉普拉斯算子为非零，采用了前面的差分公式\((6)\):
\[f_{v}(x)=m^{2} \mu \Sigma_{j}\left(\frac{\vec{v}_{j}-\vec{v}_{i}}{\rho_{j}}\right) \nabla^{2} W\left(x-x_{j}\right)
\]
4丶计算流体的压强
流体的压强是流体内部的一种力，它是维持流体不可压缩的关键。目前我们只知道流体的密度值，流体的密度与压强存在着某种联系，密度大的趋于压强必然大，因此需要根据流体的密度计算其对应的压强值。
当然我们可以直接求解关于流体压强的泊松方程\(\nabla^{2} P=\rho \frac{\nabla \vec{v}}{\Delta t}\)这样求解得到的压强是非常精确的，直接求解泊松方程在基于欧拉网格的流体模拟中比较常见。但是求解大规模稀疏矩阵的泊松方程非常耗时，因此在基于拉格朗日粒子的流体模拟比较少（但也有），一种非常廉价且效果非常不错的额方法就是采用泰特的状态方程：
\[P=B\left(\left(\frac{\rho}{\rho_{0}}\right)^{\gamma}-1\right)
\]
在上式中，\(\rho _0\)是流体静止时的密度，\(\rho\)是流体粒子当前的密度，\(\gamma\)是状态方程的指数，取值为\(7\),而\(B\)是一个缩放系数因子，与声波在流体中的传播速度有关，其计算公式如下：
\[B=\frac{\rho_{0} c_{s}^{2}}{\gamma}
\]
通过这个状态方程，我们可以计算出流体的压强值。但是存在一个问题，对于的粒子密度，我们是通过邻域粒子来计算的，这意味着流体表面的粒子因为邻域粒子较少使得其计算出来的密度值ρρ低于静止时的密度\(\rho _0\)，导致通过状态方程计算出来的压强为负数，如下图所示，从而导致流体在表面上的不正常聚集现象，这种现象类似于表面张力，但它并不是物理意义上准确的，所以我们需要消除这种现象。当计算得到的压强为负时，我们赋予为零。
5丶计算流体的梯度力
前面一步我们得到了流体粒子的压强，在压强的作用下，流体粒子从高压强区域流向低压强区域，是流体保持不可压缩的性质。因此，我们需要计算作用在流体粒子上的压强梯度力项，压力梯度力计算公式为：
\[f_{p}=-m \frac{\nabla p}{\rho}
\]
利用前面的公式\((7)\),压强梯度力离散计算公式为：
\[f_{p}=-m^{2} \Sigma_{j}\left(\frac{p_{i}}{\rho_{i}^{2}}+\frac{p_{j}}{\rho_{j}^{2}}\right) \nabla W\left(x-x_{j}\right)
\]
6丶时间步进积分
前面我们计算得到每个流体粒子的作用力合力，接着需要根据这个力计算粒子的加速度，然后在加速度的作用下更新粒子的速度值，最后在速度的作用下计算粒子的位置向量。这个步骤只需简单地利用牛顿定律即可。
7丶碰撞处理（待处理）
在获取了新得速度场和位置向量之后，我们需要做碰撞检测处理，使之不发生穿透。这里的碰撞检测采用了隐式表面的符号距离场，为每个碰撞体构建一个水平集，这方面的内容比较多，但不是这里的核心点，故不赘述。目前实现的是刚体碰撞。
8丶自适应步长
后还有个非常关键的点，就是模拟的时间步长选取。为了使得流体的密度守恒（即保持不可压缩性），我们采用了泰特的状态方程，该方程引入了两个参数，分别是指数γγ和声波在流体中的传播速度，这带来了一些时间步长的限制问题。如下图所示，假设一个流体粒子从空中落到一滩流体中，这个过程会产生一些震荡波，白色的粒子是落下的粒子以及获取到了震荡波动信息的粒子。
受限于有限的光滑核半径，在每一个时间步长内，信息传播的范围最大为光滑核半径的长度。设信息的传播速度为\(c\)，那么我们能取得最大时间步长就为\(h/c\)。在我们得物理模拟中，这个信息传播速度实际上就是声波在流体中的传播速度。为此，研究者们提出自适应的时间步长，根据当前的流体状态计算最大的时间步长，如果取超过这个最大时间步长限制，那么将会导致流体崩溃，产生不稳定的问题。时间步长的上界计算公式为：
\[\begin{aligned}
&\Delta t_{v}=\frac{\lambda_{v} h}{c_{s}} \\
&\Delta t_{f}=\lambda_{f} \sqrt{\frac{m h}{F_{\max }}} \\
&\Delta t \leq \min \left(\Delta t_{v}, \Delta t_{f}\right)
\end{aligned}
\]
其中，\(h\)是光滑核半径，\(m\)粒子质量，\(c_s\)是声波在流体中的传播速度，\(F_{max}\)是流体粒子当中受到的最大合理的长度值，\(\lambda_v\)和\(\lambda_f\)是一个缩放系数，分别取\(0.4\)和\(0.25\)。可以看到，时间步长的上界不仅仅取决于声波的传播速度，还取决于每一个不同的时刻流体所受的最大合力，因而每一次模拟都要重新计算下一个时间步长。
真实世界中的声波在流体中的传播速度很快，这导致通过上式计算得到的时间步长是非常小的。举个粒子，假设我们的光滑核半径为\(0.1m\),质量为\(0.01kg\)，声波传播速度为\(1482m/s\),在只有重力这个外力的作用下，根据上式计算可得\(\Delta t_v=0.00002699055331s\),而\(\Delta t_f=0.0007985957062\),因此最大的时间步长为\(0.00002699055331\),这意味著如果模拟\(60\)帧率的话，那么一秒需要计算\(618\)个子时间步长，耗费大量的时间和计算资源，即便是离线模拟，其代价也太大了。这个就是传统的SPH算法的缺点，被称为WCSPH（Weakly Compressible SPH），针对这个缺点，目前已经有不少研究者提出了不同的改进方法，其中比较优秀的算法就是PCISPH。
3.PCISPH
传统的SPH算法采用了泰特的状态方程来计算流体压强，避免了求解大型稀疏矩阵的泊松方程，但是受限于状态方程中的声波速度，我们不能取太大的刚度系数，否则将导致时间步长非常小，极大地耗费计算资源。因此，后来又有学者提出了传统SPH算法的改进——PCISPH（全称为Predictive-Corrective Incompressible SPH，译为预测-矫正的不可压缩SPH，简称PCISPH）。与WCSPH算法不同，PCISPH不再采用泰特方程来计算流体的压强值，而是采用了一种非线性的迭代方法，尽可能地计算一个使得流体密度守恒的压强值。
PCISPH总体算法如下图。可以看到，在求解流体的压强梯度力时，它是一个循环迭代的过程，迭代的终止条件为流体的当前密度与静止密度之差小于给定的阈值或者达到最大的迭代次数。
在每一次迭代的过程中，首先根据粒子的净合力计算加速度，并计算在当前加速度下的速度值，位置向量，这一过程为预测过程。紧接这计算前面预测得到的粒子群的密度值以及当前的密度值与静止密度的差距，并根据这个密度差，使得在这个压强的作用下，我们的粒子能够移动到一个位置上从而缩小这个密度差距。最后根据压强计算压强梯度力，将迭代到我们的合力当中，用以下一个的迭代。
迭代过程的关键一步就是如何根据密度差计算出所需的压强值，这设计到一些数学内容。设光滑核半径长度为\(h\),给定一点\(i\)在第\(t+1\)时间点的密度值通过下面的公式计算得到：
\[\begin{aligned}
\rho_{i}(t+1) &=m \Sigma_{j} W\left(x_{i}(t+1)-x_{j}(t+1)\right) \\
&=m \Sigma_{j} W\left(x_{i}(t)+\Delta x_{i}(t)-x_{j}(t)-\Delta x_{j}(t)\right) \\
&=m \Sigma_{j} W\left(d_{i j}(t)+\Delta d_{i j}(t)\right)
\end{aligned}
\]
其中\(d_{i j}(t)=x_i(t)-x_j(t)\),\(\Delta d_{i j}(t)=\Delta x_i(t)- \Delta x_j(t)\)。假设\(\Delta d_{ij}\)足够小，对\(W\left(d_{i j}(t)+\Delta d_{i j}(t)\right)\)做一阶泰勒展开：
\[\begin{aligned}
\rho_{i}(t+1) &=m \Sigma_{j}\left[W\left(d_{i j}(t)\right)+\nabla W\left(d_{i j}(t)\right) \cdot \Delta d_{i j}(t)\right] \\
&=m \Sigma_{j} W\left(x_{i}(t)-x_{j}(t)\right)+m \Sigma_{j} \nabla W\left(x_{i}(t)-x_{j}(t)\right) \cdot\left(\Delta x_{i}(t)-\Delta x_{j}(t)\right) \\
&=\rho_{i}(t)+\Delta \rho_{i}(t)
\end{aligned}
\]
上式中，\(\Delta \rho _i (t)\)是未知项，我们对其做一些展开：
\[\begin{aligned}
\Delta \rho_{i}(t) &=m \Sigma_{j} \nabla W_{i j} \cdot\left(\Delta x_{i}(t)-\Delta x_{j}(t)\right) \\
&=m\left(\Sigma_{j} \nabla W_{i j} \Delta x_{i}(t)-\Sigma_{j} \nabla W_{i j} \Delta x_{j}(t)\right) \\
&=m\left(\Delta x_{i}(t) \Sigma_{j} \nabla W_{i j}-\Sigma_{j} \nabla W_{i j} \Delta x_{j}(t)\right)
\end{aligned}
\]
令\(\Delta x_i\)w为在仅仅考虑流体压强梯度力的情况下的位移向量，则有：
\[\Delta x_{i}=\Delta t^{2} \frac{F_{i}^{p}}{m}
\]
在理想的情况下，我们的流体应该处于一种这样的平衡状态：流体粒子的压强值军等于\(\bar{p}_{i}\),流体粒子的密度值均为静止时的密度\(\rho _0\),即有：
\[F_{i}^{p}=-m^{2} \Sigma_{j}\left(\frac{\bar{p}_{i}}{\rho_{0}^{2}}+\frac{\bar{p}_{i}}{\rho_{0}^{2}}\right) \nabla W_{i j}=-m^{2} \frac{2 \bar{p}_{i}}{\rho_{0}^{2}} \Sigma_{j} \nabla W_{i j}
\]
代入得：
\[\Delta x_{i}=-\Delta t^{2} m \frac{2 \bar{p}_{i}}{\rho_{0}^{2}} \Sigma_{j} \nabla W_{i j}
\]
记在粒子\(i\)的压强作用下，邻域粒子\(j\)的位移为\(\Delta x_{j|x}\),因力是相互作用、相互对称的，则粒子\(j\)受到来自\(i\)的压强梯度力为：
\[F_{j \mid i}^{p}=m^{2}\left(\frac{\bar{p}_{i}}{\rho_{0}^{2}}+\frac{\bar{p}_{i}}{\rho_{0}^{2}}\right) \nabla W_{i j}=m^{2} \frac{2 \bar{p}_{i}}{\rho_{0}^{2}} \nabla W_{i j}
\]
在上述的压强梯度力作用下：
\[\Delta x_{j \mid i}=\Delta t^{2} m \frac{2 \bar{p}_{i}}{\rho_{0}^{2}} \nabla W_{i j}
\]
代回\(\Delta \rho _i (t)\)公式得：
\[\begin{aligned}
\Delta \rho_{i}(t) &=m\left(-\Delta t^{2} m \frac{2 \bar{p}_{i}}{\rho_{0}^{2}} \Sigma_{j} \nabla W_{i j} \cdot \Sigma_{j} \nabla W_{i j}-\Sigma_{j}\left(\nabla W_{i j} \cdot \Delta t^{2} m \frac{2 \bar{p}_{i}}{\rho_{0}^{2}} \nabla W_{i j}\right)\right) \\
&=\Delta t^{2} m^{2} \frac{2 \bar{p}_{i}}{\rho_{0}^{2}}\left(-\Sigma_{j} \nabla W_{i j} \cdot \Sigma_{j} \nabla W_{i j}-\Sigma_{j}\left(\nabla W_{i j} \cdot \nabla W_{i j}\right)\right)
\end{aligned}
\]
对上式进行变化，可得关于\(\bar{p}_{i}\)得一个表达式：
\[\begin{gathered}
\bar{p}_{i}=\frac{\Delta \rho_{i}(t)}{\beta\left(-\Sigma_{j} \nabla W_{i j} \cdot \Sigma_{j} \nabla W_{i j}-\Sigma_{j}\left(\nabla W_{i j} \cdot \nabla W_{i j}\right)\right)} \\
\beta=\Delta t^{2} m^{2} \frac{2}{\rho_{0}^{2}}
\end{gathered}
\]
而\(\Delta \rho_{i}(t)=-\rho_{\mathrm{erri}}^{t}=-\left(\rho_{i}^{t}-\rho_{0}\right)\),即当前粒子得密度值与静止密度得差。上式中的系数可以提出来提前技术，一方面是为了性能考虑，另一方面是因为表面粒子的邻域粒子数量少，其计算的系数是错的。
\[\delta=\frac{-1}{\beta\left(-\Sigma_{j} \nabla W_{i j} \cdot \Sigma_{j} \nabla W_{i j}-\Sigma_{j}\left(\nabla W_{i j} \cdot \nabla W_{i j}\right)\right)}
\]
然后\(\bar{p}_{i}=\delta \rho_{e r r i}^{t}\)，在每一次迭代中更新粒子的压强，\(p_{i}+=\bar{p}_{i}\)。PCISPH算法与WCSPH算法的不同就在计算压强梯度力上，除了更新粒子的压强值，还要更新粒子的密度值，并根据粒子的预测位置计算压强梯度力，整个迭代的流程不难理解。迭代结束之后，需要将迭代得到的压强梯度力加到粒子的合力上。
与传统的WCSPH相比，PCISPH看起来步骤更多，因为它需要迭代多次，但是其不再依赖泰特方程，从而支持更大的时间步长，因此速度比WCSPH快很多。
参考资料
博客
【笔记】SPH流体模拟基础
流体模拟Fluid Simulation：基于SPH的拉格朗日流体模拟
【Unreal从0到1】【第九章：粒子系统与物理引擎】9.1，基于Niagara的SPH可交互流体特效
参考
本构关系(Cibstutytuve Relation) 可以理解为因果关系，是关于两个场量变化关系的必然性理论。例如在可压缩流体中，流体受到的压力会引起它的体积压缩，即有本构关系。而在不可压缩流体中，流体的压力和体积压缩则没有本构关系。
