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

最新推荐文章于 2024-06-20 09:47:20 发布

置顶 0小龙虾0

最新推荐文章于 2024-06-20 09:47:20 发布

阅读量3.8k

点赞数 9

分类专栏：流体模拟(SPH)

本文链接：https://blog.csdn.net/qq_39300235/article/details/109600579

版权

流体模拟(SPH) 专栏收录该内容

13 篇文章 197 订阅

订阅专栏

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^2
    Point* 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);
        }
        
        //利用密度约束计算 deltaP
        for(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^2
    Point* 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-left
    b_density += _boundDensityContri(pi->predicted_position.x - m_sphWallBox.min.x);
    // x-right
    b_density += _boundDensityContri(m_sphWallBox.max.x - pi->predicted_position.x);
    // y-down
    b_density += _boundDensityContri(pi->predicted_position.y - m_sphWallBox.min.y);
    // y-up
    b_density += _boundDensityContri(m_sphWallBox.max.y - pi->predicted_position.y);
    // z-near
    b_density += _boundDensityContri(pi->predicted_position.z - m_sphWallBox.min.z);
    // z-far
    b_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的梯度C
        glm::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);
    }else
        pi->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的核心算法和代码就完成了，实验结果展示：