SPH（光滑粒子流体动力学）流体模拟实现六：Position Based Fluid（PBF）

PBF方法和前篇提到的PCISPH方法类似，都属于迭代矫正法。PCISPH是通过迭代预测压力，通过压力变化达到粒子的不可压缩性，而PBF则通过迭代更新位置，从而实现不可压缩性。这两者的时间步长的限制也不会那么严格，即可以有较长的时间步长也能达到不可压缩性。经过实验PBF所需要的迭代次数也要少于PCISPH。

为了直观了解二者差距，我们可以用图解作为例子，PCISPH：

我们计算出当前粒子密度和静态密度之间的差，根据差值求解压力变化量F，然后利用更新的压力求解矫正后的粒子位置。

PBF：

我们当粒子不满足约束条件时，即挤压。我们计算位移量 $\Delta P$ ，评更新其位置。

Position Based Dynamics（PBD）

了解PBF之前，我们需要了解其理论基础PBD。PBD方法中，动力学物体由N个顶点和M个约束构成。

顶点 $i\in [1,2,...,N]$ ，其质量为 $m_{i}$ ，位置为 $x_{i}$ ，速度为 $v_{i}$ 。

约束 $j\in[1,2,...,M]$ ，有以下5个属性：

1.约束的基数 $n_{j}$ ，指该约束影响的顶点个数。

2.约束函数 $C_{j}:R^{3n_{j}}\rightarrow R$ 。

3.约束索引 $\begin{Bmatrix} i_{1},i_{2},...,i_{n_{j}} \end{Bmatrix}$ ， $i_{k}\in[1,2,...,N]$ 。

4.刚度参数 $k_{j}\in[0,1]$ 。在流体中可以直接理解为约束强度。

5.约束种类，分为等式约束 $Cj(x_{i_{1}},x_{i_{2}},...,x_{i_{n_{j}}})=0$ 和不等式约束 $Cj(x_{i_{1}},x_{i_{2}},...,x_{i_{n_{j}}})\geq 0$ 。

值得注意的是，这里我们满足约束的情况下不需要矫正，而不满足情况才需要矫正。

PBD算法：

算法思路：

（1）～（3）：完成对所有顶点的初始化，并对质量求倒数，方便之后使用乘法而不是除法。

（4）～（17）：（5）：对每个粒子根据外力求解速度，这里的外力指不能转换为位置约束的力(重力)。

（6）：为了根据需要降低速度，因此引入阻尼来衰减速度。

（7）：所有的粒子利用当前的速度计算预测位置 $P_{i}$

（8）：对所有的粒子生成碰撞约束，如边界碰撞或者其他材料碰撞。

（9）～（11）：约束投影，即迭代更新预测位置 $P_{i}$ ，使其满足约束。

（12）～（15）：利用满足约束的预测 $P_{i}$ 重新求解速度，然后更新最终位置 $X_{i}$ 。

（16）：根据需要控制速度，和（6）类似，这里引入摩擦和恢复系数修改速度。

约束投影：

根据约束投影一组点意味着移动这些点以使其满足约束。且最重要的问题是线性动量和角动量的守恒。令 $\Delta P_{i}$ 为顶点i的投影位移，如果满足线性动量守恒，即：

直观理解为保持重心。

如果又满足角动量守恒，即：

$r_{i}$ 是 $P_{i}$ 到任意公共旋转中心的距离。

如果投影违反了这些约束之一，它将引入所谓的ghost force，其作用类似于外部力拖动和旋转对象。PBD则正好满足线性动量守恒和角动量守恒。PBD的 $\Delta P$ 采用约束的梯度方向进行位移。在梯度方向位移明显的保证了不会加入外部拖动和旋转力，保证了角动量守恒。且所有的约束点都进行位移后，质心也不会改变，即满足线性动量。

我们设一个约束的基数为n，该基数关联的点为 $p_{1},p_{2},...p_{n}$ ，约束函数为C，刚度系数为k。这些点位置构成的矩阵为 $p=[p_{1}^{T},p_{2}^{T},...,p_{n}^{T}]$ ，则该约束为：

我们需要求出的位移量 $\Delta P$ 能满足：

这里利用了一阶泰勒展开。我们根据上述推论，令：

将（4）带入（3）得：

这正是典型的Newton-Raphson迭代公式。然而对于粒子i来说，约束投影后对应的位移向量为：

其中：

分母是这样计算得来的：由于我们的约束是 $C_{j}:R^{3n_{j}}\rightarrow R$ ，我们的自变量可以理解为 $3*n_{j}$ 维的向量，因此该向量模的平方为。同样，我们可以看到该约束下所有点的s值是相同的。

如果我们的粒子质量是不同的呢？我们只需要考虑每个粒子i的质量 $m_{i}$ ，且其质量倒数为 $w_{i}=1/m_{i}$ ，则公式（4）会变成，计算得：

公式（6）变为：

约束例子：

我们举一个例子，如下图：

这里的距离约束为，我们对 $P_{1}$ ， $P_{2}$ 分别求梯度：

$C(P_{1},P_{2})=|P_{1}-P_{2}|-d$

$|P_{1}-P_{2}|^{2}=(P_{1}-P_{2})\cdot (P_{1}-P_{2})$

$\bigtriangledown _{P_{1}}|P_{1}-P_{2}|^{2}=2|P_{1}-P_{2}|\cdot \bigtriangledown _{P_{1}}|P_{1}-P_{2}|=\bigtriangledown _{P_{1}}((P_{1}-P_{2})\cdot (P_{1}-P_{2}))$

根据向量求导法则：，因此：

$\bigtriangledown _{P_{1}}((P_{1}-P_{2})\cdot (P_{1}-P_{2}))=\bigtriangledown _{P_{1}}(P_{1}-P_{2})\cdot (P_{1}-P_{2})+ (P_{1}-P_{2})\cdot\bigtriangledown _{P_{1}}(P_{1}-P_{2})$

$=2\cdot (P_{1}-P_{2})$

带入得：

$\bigtriangledown _{P_{1}}|P_{1}-P_{2}|=\frac{2\cdot (P_{1}-P_{2})}{2\cdot |P_{1}-P_{2}|}$

$\bigtriangledown _{P_{1}}C(P_{1},P_{2})=\frac{P_{1}-P_{2}}{|P_{1}-P_{2}|}$

同理：

$\bigtriangledown _{P_{2}}C(P_{1},P_{2})=-\frac{P_{1}-P_{2}}{|P_{1}-P_{2}|}$

带入公式（8）可得：

$s=\frac{|P_{1}-P_{2}|-d}{w_{1}+w_{2}}$

将s带入公式（9）得：

Position Based Fluids（PBF）

有了PBD的理论基础，PBF就比较简单了。就是将PBD与SPH进行结合。在不可压缩流体中，我们希望流体密度是恒定的，因此我们施加了一个密度约束：

其中 $\rho _{i}$ 是当前粒子密度， $\rho _{0}$ 是静态粒子密度。其中 $\rho _{i}$ 的计算由SPH的光滑核函数计算：

我们改文章只考虑所有粒子质量相同的情况，因此后续的公式将隐藏质量。 PBD旨在找到满足该约束的粒子位置校正∆p：

根据PBD的思想，我们的位移是沿着约束函数C的梯度方向，因此有：

我们将SPH的光滑核函数带入约束函数C，并对其求梯度得：

根据k为自身粒子还是领居粒子，这里需要分为两种情况：

在PBD中，我们知道 $\lambda$ 的求解公式为：

该值对于约束中的所有粒子都是相同的。由于约束函数C是非线性的，并且在平滑核边界处梯度逐渐消失，因此，当粒子接近分离时，该式子中的分母会导致不稳定。在PCISPH中是通过预计算的方式，它保证了在计算时领域有充足的粒子。或者，可以使用约束力混合（CFM）来规范约束。 CFM的思想是通过将一些约束力混回到约束函数中来软化约束，该方法在分母引入了 $\epsilon$ ：

这极大的加强了粒子的稳定性。

$\Delta P_{i}$ 还包括相邻粒子密度约束的校正影响，因此：

这里变号的原因是相邻粒子对当前粒子的梯度互为相反值。

XSPH人工粘性

我们这里还可以引入人工粘性（Artificial Viscosity），将流体动能转换为热能，增加数值稳定性：

其中c=0.01。

算法实现：

算法思路：

（1）～（4）：利用外力（不影响约束的力）计算速度，然后计算预测位置。

（5）～（7）：利用预测位置确定领居结构。

（8）～（19）：进入迭代矫正。（9）～（11）：对每个粒子求解 $\lambda_{i}$ 。

（12）～（15）：利用每个粒子的 $\lambda_{i}$ 求解 $\Delta P_{i}$ ，并进行碰撞检测和响应。

（16）～（18）：使用 $\Delta P_{i}$ 更新预测位置。

（20）～（24）：利用预测位置和原位置差更新速度，应用XSPH粘性，最后用预测位置更新旧位置。

代码实现

相比理论，PBF的程序实现要简单的多，我们的tick()函数如下：

    void FluidSystem::tick(bool suf){//计算外力（只包含重力），更新速度，并计算预测位置_computerExternalForces();//计算领居结构m_gridContainer.insertParticles(&m_pointBuffer);//每帧刷新粒子位置m_neighborTable.reset(m_pointBuffer.size());//重置邻接表for(int i=0;i<m_pointBuffer.size();i++){Point* pi=m_pointBuffer.get(i);//装载领居pi->kernel_self=_computeNeighbor(i);m_neighborTable.point_commit();}//迭代求解_constraintProjection();//更新位置_updatePosAndVel();//用于重建表面CalImplicitFieldDevice(m_rexSize, tem, glm::vec3((1.0/m_gridScale)/m_gridContainer.getDelta()), m_field);clearSuf();//清空表面数据m_mcMesh.CreateMeshV(m_field, tem, (1.0/m_gridScale)/m_gridContainer.getDelta(), m_rexSize, m_thre, m_vrts, m_nrms, m_face);culAll();}

根据算法思路，我们首先利用外力更新速度，并计算预测位置，_computerExternalForces()函数代码如下：

    void FluidSystem::_computerExternalForces(){for(int i=0;i<m_pointBuffer.size();i++){Point* pi=m_pointBuffer.get(i);pi->forces=glm::vec3(0.0);//外力(仅重力)pi->forces+=m_gravityDir*m_pointMass;//计算速度pi->velocity+=m_deltaTime*(pi->forces/m_pointMass)*m_unitScale;//预测位置pi->predicted_position+=pi->velocity*m_deltaTime/m_unitScale;}}

然后我们计算领居结构，值得注意的是，这里我们必须用预测的位置去计算领居结构，_computeNeighbor()函数代码如下：

float FluidSystem::_computeNeighbor(int i){float h2=m_smoothRadius*m_smoothRadius;//h^2Point* pi=m_pointBuffer.get(i);float sum=0.0;m_neighborTable.point_prepare(i);int gridCell[8];m_gridContainer.findCells(pi->predicted_position, m_smoothRadius/m_unitScale, gridCell);for(int cell=0;cell<8;cell++){if (gridCell[cell]==-1) continue;int pndx=m_gridContainer.getGridData(gridCell[cell]);while (pndx!=-1) {Point* pj=m_pointBuffer.get(pndx);if(pj==pi)sum+=pow(h2,3.0);else{glm::vec3 pi_pj=(pi->predicted_position-pj->predicted_position)*m_unitScale;float r2=pi_pj.x*pi_pj.x+pi_pj.y*pi_pj.y+pi_pj.z*pi_pj.z;if(h2>r2){float h2_r2=h2-r2;sum+=pow(h2_r2, 3.0);if(!m_neighborTable.point_add_neighbor(pndx, glm::sqrt(r2)))return sum*m_kernelPoly6;}}pndx=pj->next;}}return sum*m_kernelPoly6;
}

之后，我们正式进入迭代求解，这里我就直接使用的雅可比迭代，迭代函数_predictionCorrection()代码如下：

void FluidSystem::_constraintProjection(){int iteration=0;while (iteration<maxLoops) {//计算粒子密度和lambda，(都是使用预测位置)for(int i=0;i<m_pointBuffer.size();i++){_culDensity(i);_calculateLambda(i);}//利用密度约束计算 deltaPfor(int i=0;i<m_pointBuffer.size();i++){_calculateDeltaPi(i);_collisionHandling(i);}//更新预测位置for (int i=0; i<m_pointBuffer.size(); i++) {Point* pi=m_pointBuffer.get(i);pi->predicted_position+=pi->deltaX;}iteration++;}
}

在迭代中，我们首先需要计算粒子的 $\rho _{i}$ 和 $\lambda_{i}$ ，_culDensity()函数如下：

void FluidSystem::_culDensity(int i){float h2=m_smoothRadius*m_smoothRadius;//h^2Point* pi=m_pointBuffer.get(i);int neighborCounts=m_neighborTable.getNeighborCounts(i);float sum=pow(h2,3.0);for (unsigned int j=0; j<neighborCounts; j++) {unsigned int neighborIndex;float r;m_neighborTable.getNeighborInfo(i, j, neighborIndex, r);float r2=r*r;float h2_r2=h2-r2;sum+=pow(h2_r2, 3.0);}pi->density=sum*m_kernelPoly6*m_pointMass;//边界处理float b_density = 0.0f;// x-leftb_density += _boundDensityContri(pi->predicted_position.x - m_sphWallBox.min.x);// x-rightb_density += _boundDensityContri(m_sphWallBox.max.x - pi->predicted_position.x);// y-downb_density += _boundDensityContri(pi->predicted_position.y - m_sphWallBox.min.y);// y-upb_density += _boundDensityContri(m_sphWallBox.max.y - pi->predicted_position.y);// z-nearb_density += _boundDensityContri(pi->predicted_position.z - m_sphWallBox.min.z);// z-farb_density += _boundDensityContri(m_sphWallBox.max.z - pi->predicted_position.z);//m_kBoundsDensity=5;pi->density+= b_density*m_kBoundsDensity;
}float FluidSystem::_boundDensityContri(float dx_){dx_*=m_unitScale;static const float PI=3.1415926;if (dx_ > m_smoothRadius) {return 0.0f;}if (dx_ <= 0.0f) {return (2 * PI / 3);}return (2 * PI / 3) * pow(m_smoothRadius - dx_, 2) * (m_smoothRadius + dx_);
}

我们还考虑了边界对密度的影响。

_calculateLambda()函数代码如下：

void FluidSystem::_calculateLambda(int i){const float eps=1.0e-6;Point* pi=m_pointBuffer.get(i);const float C=max(pi->density/m_restDensity-1.0f,0.0f);if(C!=0.0){//梯度C的平方float sum_grad_C2=0.0;//pi的梯度Cglm::vec3 gradC_i=glm::vec3(0.0);int neighborCounts=m_neighborTable.getNeighborCounts(i);for(int j=0;j<neighborCounts;j++){unsigned int neighborIndex;float r;m_neighborTable.getNeighborInfo(i, j, neighborIndex, r);Point* pj=m_pointBuffer.get(neighborIndex);float h_r=m_smoothRadius-r;glm::vec3 r_V=(pi->predicted_position-pj->predicted_position)*m_unitScale;glm::vec3 grad_W=r_V*m_kernelSpiky*h_r*h_r/r;glm::vec3 gradC_j=-m_pointMass/m_restDensity*grad_W;sum_grad_C2+=glm::dot(gradC_j,gradC_j);gradC_i-=gradC_j;}sum_grad_C2+=glm::dot(gradC_i,gradC_i);//计算lambda（使用软约束，即带入eps项）pi->lambda=-C/(sum_grad_C2+eps);}elsepi->lambda=0.0f;
}

之后我们利用密度约束计算 $\Delta P$ ，并进行碰撞检测，_calculateDeltaPi()函数和_collisionHandling()函数代码如下：

void FluidSystem::_calculateDeltaPi(int i){Point* pi=m_pointBuffer.get(i);pi->deltaX=glm::vec3(0.0);int neighborCounts=m_neighborTable.getNeighborCounts(i);for(int j=0;j<neighborCounts;j++){unsigned int neighborIndex;float r;m_neighborTable.getNeighborInfo(i, j, neighborIndex, r);Point* pj=m_pointBuffer.get(neighborIndex);float h_r=m_smoothRadius-r;glm::vec3 r_V=(pi->predicted_position-pj->predicted_position)*m_unitScale;glm::vec3 grad_W=r_V*m_kernelSpiky*h_r*h_r/r;glm::vec3 gradC_j=-m_pointMass/m_restDensity*grad_W;pi->deltaX -= (pi->lambda + pj->lambda) * gradC_j;}
}void FluidSystem::_collisionHandling(int i){const float damping = 0.1;Point* pi=m_pointBuffer.get(i);for(int i=0;i<3;i++){if(pi->predicted_position[i]<m_sphWallBox.min[i]){pi->predicted_position[i]=m_sphWallBox.min[i];pi->velocity*=damping;}if(pi->predicted_position[i]>m_sphWallBox.max[i]){pi->predicted_position[i]=m_sphWallBox.max[i];pi->velocity*=damping;}}
}

这里我们迭代步数设为5次，就能达到较好的效果。完成迭代后，我们需要更新粒子的速度和位置，_updatePosAndVel()函数代码如下：

void FluidSystem::_updatePosAndVel(){for(int i=0;i<m_pointBuffer.size();i++){Point* pi=m_pointBuffer.get(i);const float invT=1.0/m_deltaTime;pi->velocity = invT*(pi->predicted_position-pi->pos)*m_unitScale;//计算人工黏力_computeXSPHViscosity(i);pi->pos = pi->predicted_position;  // p(t+1) = p(t) + v(t+1/2) dt//弹入位置数据posData[3*i]=pi->pos.x;posData[3*i+1]=pi->pos.y;posData[3*i+2]=pi->pos.z;}
}

这里我们利用预测位置计算完速度后，还需要用XSPH人工粘性矫正速度，_computeXSPHViscosity()函数代码如下：

void FluidSystem::_computeXSPHViscosity(int i){const float h2=m_smoothRadius*m_smoothRadius;Point* pi=m_pointBuffer.get(i);float sum=0.0;int neighborCounts=m_neighborTable.getNeighborCounts(i);for(int j=0;j<neighborCounts;j++){unsigned int neighborIndex;float r;m_neighborTable.getNeighborInfo(i, j, neighborIndex, r);Point* pj=m_pointBuffer.get(neighborIndex);float h_r=m_smoothRadius-r;glm::vec3 vij=pi->velocity-pj->velocity;float r2=r*r;float h2_r2=h2-r2;sum+=pow(h2_r2, 3.0);pi->velocity-=m_viscosity*(m_pointMass/m_restDensity)*vij*sum*m_kernelPoly6;}
}

到这里，我们PBF的核心算法和代码就完成了，实验结果展示：

PCISPH的稳定时间步长需要控制在 $\Delta T=0.0016$ 以内，而PBF只需控制在 $\Delta T=0.016$ 以内。而且PBF的迭代次数都是固定在2-4内就有很好的效果了，PCISPH需要迭代至误差小于阈值。