流体模拟：Position Based Fluid
作者：互联网
流体模拟：Position Based Fluid
之前博客中SPH基础中介绍的WCSPH或者PCISPH介绍的方法都是属于基于力的方法。这类方法通过计算内力（比如流体内部的粘性力和压力）和外力（比如重力和碰撞力）的合力，然后根据牛顿第二定律计算出加速度，再根据数值积分求出速度和位置信息。
基于位置动力学（Position Based Dynamics）的方法将这些物理运动通过约束表达出来，这样只需要一个求解器即可，更加方便地进行物理模拟。下面我们来讲一下PBD。
1. 位置动力学（PBD）
我们来看以下基于力和基于位置动力学的物理碰撞更新过程的对比：
可以看到基于力的碰撞检测首先在穿透发生时更新物体的速度，然后更新物体的位置，会有明显的反应延迟并且需要调整刚度（stiffness）参数（比较难调）。而基于位置动力学的碰撞检测首先只检测是否发生穿透，然后移动位置使之不发生穿透，最后再据此更新物体的速度信息。
1.1 算法过程
在PBD算法中，运动的物理由\(N\)个顶点和\(M\)个约束组成。顶点\(i\in [1,...,N]\)的质量为\(m_i\),位置为\(x_i\),速度为\(v_i\),每个约束\(j\in [1,...,M]\)有如下五个性质
- 约束的基数为\(n_j\),即第\(j\)个约束所影响的顶点数目为\(n_j\)个
- 约束函数\(C_j\):\(R^{3nj} \rightarrow R\)
- 受约束影响的顶点索引值集合\({i_1,...,i_{nj}},i_k\in [1,...,N]\)
- 每个约束都有对应的刚性参数\(k_j\in [0,1]\),这里我们可以理解为约束的强度；
- 约束分为两种，一类是等式约束即\(C_j(x_{i1},x_{i2},...,x_{in_j})=0\)，另一类是不等式约束$C_j(x_{i1},x_{i2},...,x_{in_j})\ge 0 $
给定时间步长\(\Delta t\)，PBD的运动物体模拟的算法伪代码如下所示:
在上面的算法第1步到第3步中，我们首先对顶点的位置、速度和质量倒数进行初始化，其中质量的倒数\(w_i=1/m_i\),除了可以避免冗余的除法操作外，还可以使用于静态的物体，对于静态的物体我们设为\(w_i=0\)(可以理解为质量无穷大),这样在后续的更新中都不会产生位置和速度的变化量。
第5步中的\(f_{ext}\)代表不能转换为约束形式的力（如重力），我们根据\(f_{ext}\)进行一次数值计算预测在\(f_{ext}\)的作用下的速度\(v_i\)。紧接着在第6步中我们添加阻尼的作用，阻尼可以理解为物体在运动中发生了能量耗散，从而导致速度有所衰减。
第8行主要是生成碰撞约束，物体会与周围的环境发生碰撞，例如布料落在地板上，水碰上一面墙等，这些碰撞约束在每个时间步长都发生改变，所以每一次都需要重新生成碰撞约束。有了内部约束（如不可压缩流体的密度约束）和外部约束（如流体与地面的碰撞约束）之后，我们需要根据这些约束做一个迭代求解，也就是上面伪代码中的第9行到第11行，这里我们称为约束投影步骤。从约束投影步骤我们得到服从给定约束的粒子位置，然后再第12行到第15行更新顶点粒子的速度和位置信息。最后在第16行根据摩擦系数（friction）和恢复系数（restitution）更新速度，如下图2所示。这样，一个完整的PBD物理模拟步骤就完成了。
1.2 约束投影步骤
接下来我们就针对约束投影步骤详细展开相关的内容，约束投影是PBD中的最难理解的核心部分，涉及的数学内容比较多一点。设有一个基数为\(n\)的约束，关联的粒子点为\(p_i,...,p_n\),约束函数记为\(C\),刚度系数（stiffness）为\(k\),记\(p=\left[p_{1}^{T}, \ldots, p_{n}^{T}\right]^{T}\),则等式约束函数表示为：
\[C(p)=0 \tag{1} \label{1} \]我们希望计算一个位移偏移量\(\Delta p\),使得粒子顶点在\(p+\Delta p\)出约束条件依然满足，即：
\[C(p+\Delta p)=0 \tag{2} \label{2} \]对约束条件\(C\)做一阶泰勒展开，则可得：
\[C(p+\Delta p) \approx C(p)+\nabla_{p} C(p) \cdot \Delta p=0 \tag{3} \label{3} \]为了使粒子在\(p+\Delta p\)出依然满足约束条件，我们要求方程\((3)\)得到\(\Delta p\)。PBD算法一个绝妙之处在于他将\(\Delta p\)的方向限制在约束函数的梯度方向\(\nabla C(p)\)上。如下图所示，约束\(C\)所涉及到了粒子位置回形成一个高维空间，下图为该空间中满足不同约束条件的粒子位置形成的二维等值线示意图，其中满足\(C\)约束条件的是黑色等值线。故当粒子处于下图的黑色点的位置时，不满足约束条件，如果我们沿着点所在的等值线（灰色曲线）移动，此时刚体模态（Rigid body modes）的方向与该等值线相同，新得到的位置仍然在该灰色等值线上，依然不在黑色曲线\(C\)上，即不满足约束条件。这可以理解为，约束中存在的误差依然没有得到修正。以两个粒子形成的距离约束为例，就好比同时移动了两个粒子或者该约束绕自身旋转，但是存在的误差并没有得到更正。而且这样一来还会引入系统中不存在的一种外力，导致系统动量不守恒。所以，我们希望该点的位移方向与刚体模态方向垂直，从而保证系统动量守恒，即从黑点指向红点的方向\(\nabla C\)
因此，令位移向量\(\nabla p\)为约束函数的梯度向量\(\nabla _pC\)再乘上一个标量缩放系数\(\lambda\)：
\[\Delta p=\lambda \nabla_{p} C(p) \tag{4} \label{4} \]其中的标量缩放系数\(\lambda\)，我们称之为拉格朗日乘子（Lagrange multiplier）。联立公式\((3)\)和\((4)\)我们可得：
\[\lambda=-\frac{C(p)}{\left|\nabla_{p} C(p)\right|^{2}} \tag{5} \label{5} \]再将$\lambda \(代回公式\)(4)\(我们可得\)\Delta p$的表达式：
\[\Delta p=\lambda \nabla_{p} C(p)=-\frac{C(p)}{\left|\nabla_{p} C(p)\right|^{2}} \nabla_{p} C(p) \tag{6} \label{6} \]具体到粒子\(i\)，约束投影后其对应的位移向量为：
\[\Delta p_{i}=-s \nabla_{p_{i}} C\left(p_{1}, \ldots, p_{n}\right) \tag{7} \label{7} \]其中的\(s\)如下所示，\(s\)的值对于约束函数\(C\)作用范围内地额所有点都一样：
\[s=\frac{C\left(p_{1}, \ldots, p_{n}\right)}{\sum_{j}\left|\nabla_{p_{j}} C\left(p_{1}, \ldots, p_{n}\right)\right|^{2}} \tag{8} \label{8} \]前面我们假定所有的粒子质量都相同，现在考虑粒子质量不同的情况。记粒子\(i\)的质量为\(m_i\),其质量的倒数为\(w_i=1/m_i\),则公式\((4)\)变为：
\[\Delta p_{i}=\lambda w_{i} \nabla_{p_{i}} C(p) \tag{9} \label{9} \]公式\((7)\)和公式\((8)\)变为：
\[\Delta p_{i}=-s w_{i} \nabla_{p_{i}} C\left(p_{1}, \ldots, p_{n}\right) \tag{10} \label{10} \]\[s=\frac{C\left(p_{1}, \ldots, p_{n}\right)}{\Sigma_{j} w_{j}\left|\nabla_{p_{j}} C\left(p_{1}, \ldots, p_{n}\right)\right|^{2}} \tag{11} \label{11} \]为了便于理解，接下来我们举个简单的例子应用约束投影方法。如下图所示。
上面的约束可以表示为\(C\left(p_{1}, p_{2}\right)=\left|p_{1}-p_{2}\right|-d\),位移向量记为\(\Delta p_i\)。根据约束投影方法，我们首先约束函数\(C(p_1,p_2)\)关于\(p_1\)和\(p_2\)的梯度，也就是求偏导。注意到\(C\left(p_{1}, p_{2}\right)=\left|p_{1}-p_{2}\right|-d=\left(\sqrt{\left(p_{1}-p_{2}\right)^{2}}\right)-d\),我们可以求得以下的梯度向量表达式：
\[\begin{gathered} \nabla_{p_{1}} C\left(p_{1}, p_{2}\right)=\frac{p_{1}-p_{2}}{\left|p_{1}-p_{2}\right|} \\ \nabla_{p_{2}} C\left(p_{1}, p_{2}\right)=-\frac{p_{1}-p_{2}}{\left|p_{1}-p_{2}\right|} \end{gathered} \tag{12} \label{12} \]注意，上面求到的是一个矢量，也就是我们说的梯度向量。将公式\((12)\)代入公式\((11)\)可得：
\[\begin{aligned} s &=\frac{C\left(p_{1}, \ldots, p_{n}\right)}{\Sigma_{j} w_{j}\left|\nabla_{p_{j}} C\left(p_{1}, \ldots, p_{n}\right)\right|^{2}} \\ &=\frac{\left|p_{1}-p_{2}\right|-d}{w_{1}\left|\nabla_{p_{1}} C\left(p_{1}, p_{2}\right)\right|^{2}+w_{2}\left|\nabla_{p_{2}} C\left(p_{1}, p_{2}\right)\right|^{2}} \\ &=\frac{\left|p_{1}-p_{2}\right|-d}{w_{1}+w_{2}} \end{aligned} \tag{13} \label{13} \]最后，将公式\((13)\)代入到公式\((10)\),可得约束投影计算得到的位移：
\[\begin{aligned} \Delta p_{1} &=-\frac{\left|p_{1}-p_{2}\right|-d}{w_{1}+w_{2}} w_{1} \nabla_{p_{1}} C\left(p_{1}, p_{2}\right) \\ &=-\frac{w_{1}}{w_{1}+w_{2}}\left(\left|p_{1}-p_{2}\right|-d\right) \frac{p_{1}-p_{2}}{\left|p_{1}-p_{2}\right|} \end{aligned} \]同理得\(\Delta p_2\),如下所示：
\[\Delta p_{2}=+\frac{w_{2}}{w_{1}+w_{2}}\left(\left|p_{1}-p_{2}\right|-d\right) \frac{p_{1}-p_{2}}{\left|p_{1}-p_{2}\right|} \]前面我们提到了每个约束都有对应得刚度系数\(k\),令\(k'=1-(1-5)^{1/n_s}\)去乘\(\Delta p\),这里\(n_s\)迭代之后误差为\(\Delta p\left(1-k^{\prime}\right)^{n_{s}}=\Delta p(1-k)\),与刚度系数成线性关系，而与迭代次数\(n_s\)无关。下一个时间步的位置如下所示：
\[\begin{aligned} &p_{1}^{t+1}=p_{1}^{t}+k^{\prime} \Delta p_{1} \\ &p_{2}^{t+1}=p_{2}^{t}+k^{\prime} \Delta p_{2} \end{aligned} \]1.3 约束投影求解器
前面的伪代码中我们可以看到约束投影的输入为\(M+M_{coll}\)个约束和\(N\)个点的预测位置\(p_1,...,p_N\),所需要求解的方程组是非线性非对称方程组或不等式组（碰撞约束产生的）。约束投影步骤的主要任务就是修正预测位置使新得到的校正位置满足所有约束。但是一般情况下很难找到一个适当的\(\Delta p=\left[\Delta p_{1}^{T}, \ldots, \Delta p_{n}^{T}\right]^{T}\)恰好使得所有的约束都能够同时得到满足，故我们通常采用迭代的方法按顺序依次对约束进行求解。
我们可以采用非线性高斯-赛德尔（Non-Linear Gauss-Seidel，简称NGS）迭代方法。高斯赛德尔（Gauss-Sedel，简称GS）迭代方法只能求解线性方程组，NGS在依次求解德基础上，加入了约束投影求解这一非线性操作。与雅可比迭代方法（Jacobi method）不同，NGS求解器在一次迭代中对于顶点位置的修正立即被应用到下一个约束求解中，这样的好处就是显著加快了收敛速度。
但是NGS虽然稳定且容易实现，但是该方法收敛速度依然不是很快，不宜并行化。而且取决于约束求解顺序，且不宜并行化。
2 基于位置动力学的流体模拟（PBF）
前面部分主要介绍了Position Based Dynamics算法相关的内容，接下来我们就看看如何将其PBD算法应用到流体模拟当中。
2.1 不可压缩约束和拉格朗日乘子
在不可压缩流体模拟中，我们希望粒子\(i\)的密度尽量与静态密度\(\rho_0\)相同，记\(\rho_i=\rho_0\),因此，我们需要对每一个流体粒子都施加一个常量密度约束，PBF的常量密度约束如下所示：
\[C_{i}\left(p_{1}, \ldots, p_{n}\right)=\frac{\rho_{i}}{\rho_{0}}-1 \tag{14} \label{14} \]我们记粒子\(i\)的位置为\(p_i\),\(p_i,...,p_n\)使与粒子\(i\)相邻的粒子。我们可以看到当密度约束\(=0\)时，有\(\rho_i=\rho_0\),此时流体的体积既不压缩也不膨胀，从而保证了流体的不可压缩条件。
流体粒子\(i\)的密度根据SPH方法的插值计算得到：
\[\rho_{i}=\sum_{j} m_{j} W\left(p_{i}-p_{j}, h\right) \tag{15} \label{15} \]在公式\((15)\)中，\(m_j\)时邻居粒子\(j\)的质量，\(h\)是指定的光滑核半径。将公式\((15)\)代入公式\((14)\)，我们有：
\[C_{i}\left(p_{1}, \ldots, p_{n}\right)=\frac{\Sigma_{j} m_{j} W\left(p_{i}-p_{j}, h\right)}{\rho_{0}}-1 \tag{16} \label{16} \]在公式\((15)\)的密度计算中，PBF方法采用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} \tag{17} \label{17} \]但是在计算密度的梯度时，却又采用了Spiky核函数：
\[W_{spkiy}(r, h)=\frac{15}{\pi h^{6}} \begin{cases}\left(1-\frac{r}{h}\right)^{3} & 0 \leq r \leq h \\ 0 & \text { otherwise }\end{cases} \tag{18} \label{18} \]对公式\((18)\)求关于\(r\)的导数（注意，\(|r|=\sqrt{r^{2}}\)，不能直接对\(|r|\)求导），从而流体粒子密度的梯度如下所示：
\[\nabla W_{\text {spiky }}(r, h)=-\frac{45}{\pi h^{6}} \begin{cases}(h-|r|)^{2} \frac{r}{|r|} & 0 \leq|r| \leq h \\ 0 & \text { otherwise }\end{cases} \tag{19} \label{19} \]因此，粒子\(i\)的约束函数\((16)\)是一个关于\(p_1,...,p_n\)的非线性方程组，\(C_{i}\left(p_{1}, \ldots, p_{n}\right)\)=0,所有粒子的约束组成了一个非线性方程组。在PBF方法中，我们只考虑粒子质量相同的情况，我们省去公式\((15)\)和\((16)\)中的质量\(m_j\),即：
\[\begin{eqnarray*} \rho_{i}=\Sigma_{j} W\left(p_{i}-p_{j}, h\right) \tag{20} \label{20} \\ C_{i}\left(p_{1}, \ldots, p_{n}\right)=\frac{\Sigma_{j} W\left(p_{i}-p_{j}, h\right)}{\rho_{0}}-1 \tag{21} \label{21} \end{eqnarray*} \]然后求约束函数\(C_i\)关于\(p_k\)的梯度如下,其中\(k\in {1,2,...,n}\):
\[\nabla_{p_{k}} C_{i}=\frac{1}{\rho_{0}} \Sigma_{j} \nabla_{p_{k}} W\left(p_{i}-p_{j}, h\right) \tag{22} \label{22} \]显然，针对\(k\)的不同，分为两种情况。当\(k=i\)也就是粒子本身时候，连加符号中\(W\)均为关于\(p_k\)的函数；当\(k=j\)即邻居粒子的时候，只有\(W(p_i-p_k,h)\)才有意义，其他相对于\(p_k\)来说都是常量，故导数为0（注意用到了求导的链式法则）：
\[\nabla_{p k} C_{i}=\frac{1}{\rho_{0}} \begin{cases}\Sigma_{j} \nabla_{p k} W\left(p_{i}-p_{j}, h\right) & \text { if } k=i \\ -\nabla_{p k} W\left(p_{i}-p_{j}, h\right) & \text { if } k=j\end{cases} \tag{23} \label{23} \]既然求出了约束函数的梯度，我们就把它应用到前面提到的拉格朗日乘子的计算公式中，联立公式\(*(5)\)和\((23)\)，我们可得拉格朗日乘子：
\[\lambda_{i}=-\frac{C_{i}\left(p_{1}, \ldots, p_{n}\right)}{\Sigma_{k}\left|\nabla_{p k} C_{i}\right|^{2}} \tag{24} \label{24} \]2.2 软约束（Soft constraint）与混合约束力（Constraint Force Mixing, CFM）
如果一个约束条件不能被违背，我们称之为硬约束；而能一定程度上被违背的约束称为软约束。在理想的情况下，我们都希望约束始终是硬约束，但是由于误差或者数值方法的不稳定等原因，我们有时不得不向软约束妥协。
在\(PBF\)中，当\(|r|=h\),粒子\(i\)与粒子\(j\)之间的距离等于光滑核半径时，粒子\(i\)和粒子\(j\)处于即将分离的状态。注意观察公式\((19)\)的密度梯度计算公式，此时\(\nabla W_{spiky}(r,h)=0\)。若所有的邻居粒子与粒子\(i\)都处于这种状态，那么必将导致约束函数的梯度即公式\((22)\)取值为0：
\[\nabla_{p k} C_{i}=\frac{1}{\rho_{0}} \Sigma_{j} \nabla_{p k} W\left(p_{i}-p_{j}, h\right)=0 \]从而导致公式\((24)\)中的分母\(\Sigma_{k}\left|\nabla_{p k} C_{i}\right|^{2}\)为0，出现除零错误，这将导致PBF方法出现潜在的不稳定性。为了解决这个问题，PBF采用混合约束的方法，使密度硬约束转变成软约束。具体的做法就是将根据密度函数求解得到的约束力在加入到原始的约束函数中，这里的在PBF的常量密度约束中得到的拉格朗日乘子\(\lambda\)又类似的作用，故将\(\lambda\)加入到初始的约束方程（即公式\((3)\)）:
\[C(p+\Delta p) \approx C(p)+\nabla C^{T} \nabla C \lambda+\epsilon \lambda=0 \tag{25} \label{25} \]公式\((25)\)中的$\epsilon \(是松弛参数，可以由用户指定。引入公式\)(25)\(后，拉格朗日乘子的计算公式\)(24)$就变成：
\[\lambda_{i}=-\frac{C_{i}\left(p_{1}, \ldots, p_{n}\right)}{\Sigma_{k}\left|\nabla_{p k} C_{i}\right|^{2}+\epsilon} \tag{26} \label{26} \]从而可得粒子\(i\)在经过上述约束投影后对应的位移向量（包括自身密度约束以及邻居粒子密度约束共同作用的结果。注意，这里对应上面的公式\((4)\)，结合公式\((23)\)):
\[\begin{aligned} \Delta p_{i} &=\lambda_{i} \nabla_{p_{i}} C_{i}+\Sigma_{j} \lambda_{j} \nabla_{p_{j}} C_{i} \\ &=\frac{1}{\rho}_{0} \Sigma_{j} \lambda_{i} \nabla_{p_{i}} W(r, h)+\left(-\frac{1}{\Sigma_{j}} \lambda_{j} \nabla_{p_{j}} W(r, h)\right) \\ &=\frac{1}{\rho_{0}} \Sigma_{j} \lambda_{i} \nabla_{p_{i}} W(r, h)+\frac{1}{\rho_{0}} \Sigma_{j} \lambda_{j} \nabla_{p_{i}} W(r, h) \\ &=\frac{1}{\rho_{0}} \Sigma_{j}\left(\lambda_{i}+\lambda_{j}\right) \nabla_{p_{i}} W(r, h) \end{aligned} \tag{27} \label{27} \]2.3 拉伸不稳定性
PBF采用SPH的方法来计算流体粒子的密度，但是该方法通常需要30~40个邻居粒子才能使密度求值结果趋于静态密度。在邻居粒子数量较少的情况下，通过该方法计算得到的流体密度低于静态密度，由此会造成流体内部压强为负数，原本粒子间的压力变为吸引力，使得流体粒子凝聚在一起，导致流体表面的模拟效果失去真。PBF采用了一种人工排斥力的计算模型，当流体粒子距离过近时该排斥力会使它们分开，避免产生粒子聚集的现象。如下图所示：
在公式\((24)\)的基础上，加入一个排斥项（repulsice term）\(S_{corr}\):
\[\Delta p_{i}=\frac{1}{\rho_{0}} \Sigma_{j}\left(\lambda_{i}+\lambda_{j}+s_{\text {corr }}\right) \nabla_{p_{j}} W\left(p_{i}-p_{j}, h\right) \tag{28} \label{28} \]其中的\(S_{corr}\)计算方式如下：
\[s_{c o r r}=-k\left(\frac{W\left(p_{i}-p_{j}, h\right)}{W(\Delta q, h)}\right)^{n} \tag{29} \label{29} \]公式\((29)\)中，\(\Delta q\)表示到粒子\(i\)的一个固定距离，通常取\(|\Delta q|=0.1h,...,0.3h\),\(h\)即前面提到的光滑核半径。此外，公式中的\(k\)可以看作表面张力参数，取值\(k=0.1\),而\(n=4\)。公式\((28)\)中的排斥项会使得流体粒子的密度稍微低于静态密度，从而产生类似于表面张力的效果，使得流体表面的的粒子分布均匀。通过这个排斥项，我们不再需要硬性规定流体的邻居数量必须在30~40个，进一步提升算法的流体模拟效率。
2.4 涡轮控制和人工粘性
由于数值耗散，PBD的方法会引入额外的阻尼，使得整个系统的能量损耗太快，导致本来应该本来一些涡流细节迅速消失。在这里，PBF通过涡轮控制方法向整个系统重新注入能量：
\[f_{i}^{\text {vorticity }}=\epsilon\left(N \times \omega_{i}\right) \tag{30} \label{30} \]上述公式中，\(N=\frac{\eta}{|\eta|}, \eta=\nabla|\omega|_{i}\),而流体粒子的旋度\(w_i\)计算公式如下：
\[\omega_{i}=\nabla \times v=\Sigma_{j}\left(v_{j}-v_{i}\right) \times \nabla_{p_{j}} W\left(p_{i}-p_{j}, h\right) \tag{31} \label{31} \]\[\eta=\sum_{j} \frac{m_{j}}{\rho_{j}}\left\|\omega_{j}\right\| \nabla_{i} W_{i j} \]涡轮控制方法的基本思路就是：通过添加一个体积力\(f_{i}^{\text {vorticity }}\)(在算法的第一步)，在旋度粒子（可直观理解为比周围粒子旋转快的粒子，旋度\(w_i\)指向粒子\(i\)的旋转轴）处加速例子的旋转运动，通过这种方式来增加系统的旋度细节。公式\((30)\)中的\(\epsilon\)用于控制涡轮控制力的强度。
最后，PBF方法采用XSPH的粘度方法直接更新速度，从而产生粘性阻尼。人工粘性除了可以增加模拟的数值稳定性，还可以消除非物理的流体振荡。拉格朗日流体模拟方法中，人工粘性本质上会对流体粒子的相对运动产生阻尼作用，使流体的动能转化为热能：
\[v_{i}^{\text {new }}=v_{i}+c \Sigma_{j}\left(v_{i}-v_{j}\right) \cdot W\left(p_{i}-p_{j}, h\right) \tag{32} \label{32} \]在流体模拟中，我们取公式\((32)\)中的\(c=0.01\)
3. PBF算法实现
PBF算法的总体框架就是按照前面提到的PBD算法，只是经典PBD算法采用了顺序高斯-赛德尔（Sequential Gauss-Seidel，SGS）迭代求解，而SGS不容易被GPU并行化，因此基于CUDA实现的PBF求解器使用了雅克比（Jacobi）迭代方法并行求解。
参考资料
【深入浅出 Nvidia FleX】(1) Position Based Dynamics
【深入浅出 Nvidia FleX】(2) Position Based Fluids
标签：right,Based,nabla,frac,Fluid,约束,Delta,Position,left 来源： https://www.cnblogs.com/Ligo-Z/p/16263141.html