PBF流体模拟(WIP)

## 简介

> 注意：PBF目前仍然处于开发阶段，不能实用。本文档仅作为开发记录。文档将会动态更新。

PBD模块中的流体模拟目前默认称之为PBF_BVH
PBF的意思是Position Based Fluid
BVH的意思是我们目前采用BVH进行邻域搜索。

它是Miles Macklin 和Matthias Muller 2013年发表在TOG上的文章所公开的。具体请参考：

【论文】
https://mmacklin.com/pbf_sig_preprint.pdf
【Youtube展示效果视频】
https://www.youtube.com/watch?v=F5KuP6qEuew

zeno默认的FLIP求解器相比，它最大的特点就是不需要建立网格。类似houdini中的vellum grain。因为是纯粒子法，所以**输入数据只是粒子位置pos，输出数据也只有粒子位置pos**。其他所有数据都只是辅助性的。

我们先简要介绍**用法**，然后介绍**原理**。

## 用法
我们先要给出输入数据，也就是初始粒子的位置。然后连接PBF_BVH即可。
输入：prim中只需要pos。
输出：outPrim中只有pos。
参数：暂定，将会调整。
（pos数据结构用法与格式几乎与`std::vector<vec3f>`一致，并且两者可以隐式转换。）
![](/media/202210/2022-10-25_112842_306491.png)

运行效果：
由于目前参数为ad-hoc，所以场景设置有些大。需要在PBF_BVH上使用alt+F缩放到合适的大小。
![](/media/202210/PBF_BVH_1666668855.gif)

### 表面重建
参考【zeno水打柠檬】示例
http://doc.zenustech.com/project-3/doc-191/

我们为其进行表面重建。
核心节点为**SmoothLiquidSDF**。

对于从源码构建者，需要开启编译选项。
- ZENO_WITH_FastFLIP
- ZENO_WITH_zenvdb
- ZENO_WITH_ZenoFX
并且为了开启这些编译选项，请先安装openvdb。

由于SmoothLiquidSDF节点需要输入的数据格式为vdbPoint，所以我们先要把pos数据转换为vdbPoint格式。那么就先采用PrimToVDBPointDataGrid节点。该节点的Dx可以暂时给粒子直径，即6.0

然后连接SmoothLiquidSDF即可。由于该节点输出的是SDF，所以我们还需要转换回prim。也就是SDFToPrim节点。

对于SmoothLiquidSDF，还有一个必要输入数据为Dx。就是网格大小。这个值越小，重建出来的表面越细致。我们可以先给1.0

总结起来就是

![](/media/202210/2022-10-25_115333_873025.png)

表面重建后的结果

![](/media/202210/2022-10-25_115808_863409.gif)
## 原理

PBF的原理主要参考论文【Macklin2012】
https://mmacklin.com/pbf_sig_preprint.pdf

如同所有PBD方法一样，它也是分为三步
1. 外力求解（preSolve）
2. 约束求解（solveConstraint）
3. 更新速度(postSolve)

它所定义的约束是：密度约束。

如Macklin2012公式1所示

![](/media/202210/2022-10-27_183119_690828.png)
其中的p就是各个粒子的位置。（1,2,...n实际上是一个邻域内的粒子，而非所有的粒子。）

即当前密度除以静息密度 - 1.0

实际上可以将其看成就是密度误差。

因此，PBF的核心就是保证**密度的恒定性**，密度的恒定性也就是不可压缩性。（从直观上讲，粒子既不过密也不过稀疏）

为了计算C，首先应计算密度rho_i。只需要借助SPH插值公式即可。W是核函数，充当的就是插值函数的角色。m是粒子质量，如无特别需要可以统一设定为一个标量。

![](/media/202210/2022-10-27_183511_331824.png)

那么，我们还需要密度的梯度，从而在每次约束投影的时候，修正位置以保证密度的守恒性。

密度的梯度也借助SPH插值公式。SPH插值的特点是可以将物理量上的gradient转移到核函数W上。而核函数的梯度不过是个根据粒子距离计算的helper function罢了。Macklin2012中取W核函数为poly6，而gradW核函数为spikyGradient。

![](/media/202210/2022-10-27_183749_372818.png)

这里的gradC是需要分类的。也就是k为i本身，和k为邻居j。这一点和SPH的pressure gradient force的计算方法是一致的，也是为了保证对称性。

有了C和gradC，我们就可以计算PBD框架中的拉格朗日乘数lambda了。

![](/media/202210/2022-10-27_184238_564866.png)

其中epsilon是为了防止分母为0的很小的数字，可以取1e-6.

最后，仍然按照PBD框架模板，我们只要根据lambda和gradW修正位置p即可实现约束投影。

![](/media/202210/2022-10-27_184453_702016.png)

上面的公式计算出了为了保证密度恒定性应该修正的粒子位置，也就是约束投影。其中sCorr是为了防止particle deficiency造成的lumping现象而增加的人工压力。如果不加，也可以模拟，但会造成artifacts。

这个人工压力参数由下面的公式计算
![](/media/202210/2022-10-27_184629_696966.png)
Macklin2012中推荐的参数组合为：n=4, delta_q=0.3h, k=0.1

还可以使用其他方法防止lumping，比如clamping密度到rho0以上。也就是不允许rho_i小于rho0，否则就强行设置为rho0.