泛动引擎 > 技术专题 > 粒子系统 > 投影法

投影法

1. 简介

投影法是通过求解流体压强泊松方程来获取压强，来消除速度散度误差或密度误差的方法。该方法相比于PBF法虽然计算效率较低，需要使用更大的邻域半径以及更多的迭代次数才能保证仿真的稳定性，但该方法具有更高精度和仿真质量。

2. 方法

当流体被压缩时，流体个别位置会产生速度散度误差或密度误差，此时流体需要产生相应的压强来抵抗压缩，消除散度误差和密度误差以保证不可压缩性条件。而压强泊松方程则可以根据散度误差、密度误差来计算流体的压强。以下是压强泊松方程的推导过程：

Step 1. 忽略掉粘性作用、体积力作用的离散化后的流体Navior Stokes方程形式为： $$ \mathbf{v}_i^{new} = \mathbf{v}_i^{old} - \frac{\delta t}{\rho_0} \nabla_i p
$$ （公式1）
Step 2. 由于$\mathbf{v}_i^{old}$是不满足不可压缩性条件的速度场，即：$\nabla\mathbf{v}_i^{old}\ne0$ ,而待求的满足不可压缩性条件的速度场其散度应为零，即：$\nabla\mathbf{v}_i^{new}=0$。则对公式1的等式两端同时取散度，即可得到如下的压强泊松方程。 $$ \nabla_i \cdot \frac{\delta t}{\rho_0} \nabla p =-\nabla \cdot \mathbf{v}_i^{old} $$ (公式2)
Step 3. 使用粒子对压强泊松方程中的二阶微分算符（拉普拉斯算符$\nabla\cdot\nabla_i$）与一阶微分算符（散度$\nabla\cdot$）进行离散化，将压强泊松方程转化为线性系统。
Step 4. 由于压强泊松方程只考虑了速度散度误差，而没有考虑密度误差，因此需要对泊松方程的等式右端项进行修正，一种推荐的方法即是： $$ \nabla_i \cdot \frac{\delta t}{\rho_0} \nabla p =-\nabla \cdot \mathbf{v}_i^{old} + \Lambda *\max(\frac{\rho_i}{\rho_0}-1,0) $$ (公式3)

其中系数$\Lambda$为常数。
Step 5. 带入流体的边界条件（包括气液边界条件和固液边界条件），对原压强泊松方程的系数矩阵进行修正，该泊松方程会转化为对称正定的线性系统。
Step 6. 使用共轭梯度法、雅克比法等迭代算法可快速求解这一离散化后的压强泊松方程，获得当前帧中各粒子的压强值。
Step 7. 将粒子压强带入到公式1中，修正流体粒子的速度值，获得满足不可压缩性条件的粒子速度$\mathbf{v}_i^{new}$。
Step 8. 根据粒子速度$\mathbf{v}_i^{new}$来更新粒子位置，即：$\mathbf{x}_i^{new} = \mathbf{x}_i^{old}+ \mathbf{v}_i^{new}*\delta t$

以上即为基于SPH投影法的实现流程。投影法中的压强二阶微分项和速度散度项离散化方法种类繁多，不同离散化方法也有不同的特性。目前PeriDyno中选用了变分框架下使用staggered粒子的微分离散化方案来构建投影法求解器，相比其它方法，该离散化方案具有更好的收敛性、稳定性，能够有效抑制SPH的零能问题和拉伸失稳问题的影响。具体的实现方法可参考文献。

3. 实现

PeriDyno中，投影法求解器的源文件为：

"Engine/Dynamics/Module/VariationalApproximateProjection.cu";
"Engine/Dynamics/Module/VariationalApproximateProjection.h";

投影法求解器的输入为不满足不可压缩性条件的流体粒子位置和速度场，输出为满足不可压缩性条件、同时可消除密度偏移误差的的速度场。

a. 变量：

Attribute (粒子属性)
```
  DEF_ARRAY_IN(Attribute, Attribute, DeviceType::GPU, "");
```
为了准确地对流体施加固壁边界条件，该投影法中的固体需要使用Ghost粒子进行采样。而为了区别固体粒子与流体粒子，则需要引入属性数组来对粒子进行标记。
Normal (固体Ghost粒子法向量)
```
  DEF_ARRAY_IN(Coord, Normal, DeviceType::GPU, "");
```
为了准确地对流体施加滑移（Free-slip）或无滑移(No-slip)流固边界条件，固体Ghost粒子需要具有法向量参数。
NeighborList (邻域粒子列表)
```
DEF_ARRAYLIST_IN(int, NeighborIds, DeviceType::GPU, "");
```
为了实现粒子各项物理量及其微分量的插值计算，计算过程需要使用到粒子的邻域列表。该邻域列表需要在每一帧计算开始时预先由邻域查找模块生成，邻域查找模块的位置为：
```
"Engine\Framework\Topology\NeighborPointQuery.h"
"Engine\Framework\Topology\NeighborPointQuery.cu"
```

b. 函数：

泊松方程离散化相关函数
```
VC_ComputeAlpha(...)
VC_ComputeDiagonalElement(...)
VC_ComputeAx(...)
```
以上三个函数为泊松方程（公式2）中压强二阶微分项离散化的相关函数。
```
VC_ComputeDivergence(...)
```
该函数为泊松方程速度散度的离散化相关函数，用于计算粒子的速度散度误差。
```
VC_CompensateSource(...)
```
该函数用于补偿压强泊松方程的源项（即：公式2泊松方程的右端项），对应于第2节“方法”中的step 3。SPH投影法如果仅仅修正速度散度误差并不能有效地保证不可压缩性条件，因为流体的密度会逐渐发生漂移，引入严重的错误。该函数能够有效地避免流体的体积丢失问题。
```
VC_UpdateVelocity1rd(...)
```
该函数根据粒子压强来计算粒子的压强梯度，并更新流体粒子的速度（对应于第2节“方法”中的step 7 与 step 8 ），计算过程中考虑了固壁边界对于流体的作用。

参考文献

[1] He X, Wang H, Wang G, et al. A Variational Staggered Particle Framework for Incompressible Free-Surface Flows[J]. arXiv preprint arXiv:2001.09421, 2020.