SPH（光滑粒子流体动力学）流体模拟实现五：PCISPH

最新推荐文章于 2024-03-05 19:36:37 发布

0小龙虾0

最新推荐文章于 2024-03-05 19:36:37 发布

阅读量6.1k

点赞数 15

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

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

版权

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

13 篇文章 178 订阅

订阅专栏

SPH（光滑粒子流体动力学）流体模拟实现五：PCISPH

我们知道真实的液体是不可压缩的，但我们在计算机中离散的计算流体运动，在一定的时间步长内，用标准的SPH方法求解，在在粒子聚集处容易发生挤压，造成压缩。有两种常用的方法模拟不可压缩性：1.在WCSPH(弱可压缩SPH)中，利用刚性状态方程（EOS）建模压力。2.通过求解压力泊松方程实现不可压缩性。但这两种方法都有很昂贵的计算费用。

文章“Predictive-Corrective Incompressible SPH”中，提出了一种预测矫正的方法，来使粒子达到不可压缩性，其性能上相比传统两种方法，更加的高效。

PCISPH模型

SPH概述

在拉格朗日粒子描述下, 控制流体运动的偏微分方程 Navier-Stokes 方程可表示为：

SPH方法的核心思想是以离散化粒子的形式来表征连续的场, 并对场量使用积分近似的方式进行计算。位置在xi粒子i的场量：

第i个粒子的密度计算公式为：

压力场直接从Navier-Stokes方程式推导而得：

PCISPH算法

在PCISPH方法中，速度和位置会及时更新，并估计新的粒子密度。然后，对于每个粒子，计算出参考密度的预测变化，并用于更新压力值，压力值又进入压力的重新计算。此过程一直迭代到收敛，即直到所有粒子密度波动均小于用户定义的阈值η(例如1％)。在完成矫正后，我们再更新速度和位置。详细算法流程图如下：

该算法总体思路如下：

1.计算每个粒子的邻居信息，并记录在邻接表内。

2.计算出了压力之外的所有其他力(黏力，重力)。

3.执行矫正循环：执行1).，2).，3).。

1).预测所有粒子新的的速度和位置。

2).预测所有粒子新的密度，以及计算新密度和旧密度之间的差值。

3).预测所有粒子的压力。

4.完成矫正后，利用矫正的压力计算粒子速度和位置。

压强导数

我们的算法是根据预测的密度计算新的压力值，因此我们需要找到它们之间变化的关系。我们的目的是找到一个压强p，该压强以这样一种方式改变粒子的位置，即预测的密度与参考密度相对应。

对于给定的内核平滑长度h，使用SPH密度求和公式计算时间点t + 1处的密度：

其中 $d_{ij}(t)=x_{i}(t)-x_{j}(t)$ ， $\triangle d_{ij}(t)=\triangle x_{i}(t)-\triangle x_{j}(t)$ 。假定 $\triangle d_{ij}(t)$ 非常小，我们可以用一阶泰勒展开近似：

带入 $\rho _{I}(t+1)$ 得：

我们令，因此密度增量公式为：

我们假设邻居具有相等的压强 $\widetilde{p_{i}}$ 且密度为静止密度 $\rho _{0}$ (根据不可压缩性条件)，则结果是：

PCISPH算法在每次迭代校正时只矫正d流体粒子的压力，通过迭代计算出的压力让粒子之间不至于靠的过近(粒子之间距离太近可以理解为流体可压缩)。在只考虑压力 $F_{i}^{p}$ 的情况下，根据时间积分可以计算出粒子在该压力作用下的位移：

粒子i在获得粒子j的压力同时，也会对领居粒子j施加反作用力，因此：

同样，由于粒子i导致粒子j产生的位移为：

将位移增量带入密度增量公式：

因此 $\widetilde{p_{i}}$ 为：

其中 $\beta$ 为：

$\widetilde{p_{i}}$ 公式表示了我们迭代过程中，需要不断变化 $\triangle \rho _{i}(t)$ 的密度值，从而满足粒子不可压缩性。而改变 $\triangle \rho _{i}(t)$ 的密度需要改变 $\widetilde{p_{i}}$ 的压强。我们每次预测的粒子密度为： $\rho _{i}^{*}$ 。预测密度和静止密度之间的误差为： $\rho _{err}^{*}=\rho _{i}^{*}-\rho _{0}$ 。而我们的目标是让粒子密度为 $\rho _{0}$ ，这需要的密度改变量为 $-\rho _{err}^{*}$ 。因此，我们需要施加的压强值为：

这个压强计算公式在邻居粒子数目不足的时候会导致计算错误，解决办法是进行一次预计算，即在流体粒子周围充满邻居粒子的情况下计算。我们可以直接在初始化流体时计算一次系数 $\delta$ 即可：

因此，我们的压强改变量 $\widetilde{p_{i}}$ 为：

由于只要不满足不可压缩条件，我们就重复进行预测校正步骤，因此，我们需要在迭代中不断矫正压强：

算法实现

我们之前提到，在计算系数 $\delta$ 时需要在流体初始化时计算。因此我们使用函数_computeGradWValues()和_computeDensityErrorFactor()计算该系数：

    void FluidSystem::_init(unsigned int maxPointCounts, const fBox3 &wallBox, const fBox3 &initFluidBox, const glm::vec3 &gravity){        
        m_pointBuffer.reset(maxPointCounts);
        m_sphWallBox=wallBox;
        m_gravityDir=gravity;
        
        //根据粒子间距计算粒子质量
        m_pointMass=m_restDensity*pow(m_pointDistance,3.0);
        //设定流体块
        _addFluidVolume(initFluidBox, m_pointDistance/m_unitScale);

        //MarchingCube算法属性
        m_mcMesh=rxMCMesh();
        m_gridContainer.init(wallBox, m_unitScale, m_smoothRadius*2.0, 1.0,m_rexSize,m_gridScale);//设置网格尺寸(2r)
        m_field=new float[(m_rexSize[0]+1)*(m_rexSize[1]+1)*(m_rexSize[2]+1)]();
        posData=std::vector<float>(3*m_pointBuffer.size(),0);
        m_hPosW.resize(3*m_pointBuffer.size(),0.0);
        
        m_gridContainer.insertParticles(&m_pointBuffer);//每帧刷新粒子位置
        m_neighborTable.reset(m_pointBuffer.size());//重置邻接表
        
        //计算W的梯度
        _computeGradWValues();
        //计算因子
        _computeDensityErrorFactor();
    }

其中_computeGradWValues()代码如下：

void FluidSystem::_computeGradWValues(){
        float h2=m_smoothRadius*m_smoothRadius;//h^2
        const int numP=m_pointBuffer.size();
        for (int i=0; i<numP; i++) {
            Point *pi=m_pointBuffer.get(i);
            pi->sum_grad_w=glm::vec3(0.0f);
            pi->sum_grad_w_dot=0.0f;
        }
        for (int i=0; i<numP; i++) {
            Point *pi=m_pointBuffer.get(i);
            pi->kernel_self=_computeNeighbor(i);
            m_neighborTable.point_commit();
            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);
                //需要用预测的位置去重新计算距离和核值
                glm::vec3 pi_pj=(pi->pos-pj->pos)*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;
                    r=pow(r2,0.5f);
                    float h=m_smoothRadius;
                    
                    glm::vec3 gradVec=(pi->pos-pj->pos)*m_kernelSpiky/r*(h-r)*(h-r);
                    pi->sum_grad_w+=gradVec;
                    pj->sum_grad_w-=gradVec;
                    
                    pi->sum_grad_w_dot+=glm::dot(gradVec,gradVec);
                    pj->sum_grad_w_dot+=glm::dot(-1.0f*gradVec,-1.0f*gradVec);
                    
                }
            }
        }
    }

我们在其中计算所有粒子的和。

之后我们在函数_computeDensityErrorFactor()中计算系数 $\delta$ ，代码如下：

void FluidSystem::_computeDensityErrorFactor(){
        float h2=m_smoothRadius*m_smoothRadius;
        int maxNeighborIndex = 0;
        int maxNeighs = 0;
        const int numParticles = m_pointBuffer.size();
        for (int i=0; i<numParticles; i++) {
            Point *pi=m_pointBuffer.get(i);
            int neighborCounts=m_neighborTable.getNeighborCounts(i);
            int numNeighbors=0;
            //预测密度计算
            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);
                //需要用预测的位置去重新计算距离和核值
                glm::vec3 pi_pj=(pi->pos-pj->pos)*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){
                    numNeighbors++;
                }
            }
            //获取邻居最多的粒子，和邻居个数
            if (numNeighbors>maxNeighs) {
                maxNeighs=numNeighbors;
                maxNeighborIndex=i;
            }
        }
        //获取邻居最多的粒子
        Point* pmax=m_pointBuffer.get(maxNeighborIndex);
        
        //计算新的压力根据密度误差
        float restVol=m_pointMass/m_restDensity;
        float preFactor=2.0*restVol*restVol*m_deltaTime*m_deltaTime;
        float gradWTerm=glm::dot(pmax->sum_grad_w*(-1.0f),pmax->sum_grad_w)-pmax->sum_grad_w_dot;
        float divisor=preFactor*gradWTerm;
        m_densityErrorFactor=-1.0/divisor;
        const float factor = m_deltaTime / 0.0001f;
        float densityErrorFactorParameter = 0.05 * factor * factor;
        m_densityErrorFactor*=1.0*densityErrorFactorParameter;
            
    }

值得注意的是，我们还在 $\delta$ 系数之前乘上densityErrorFactorParameter影响系数，该系数受 $\Delta t$ 影响。如果不乘上该系数就会出错，论文里并没有明确提出。

计算完成预计算的 $\delta$ 系数后，我们的tick()函数如下：

    void FluidSystem::tick(bool suf){
        //求解领居结构
        m_gridContainer.insertParticles(&m_pointBuffer);//每帧刷新粒子位置
        m_neighborTable.reset(m_pointBuffer.size());//重置邻接表
        
        //计算外力
        _computerExternalForces();
        //预测矫正
        _predictionCorrection();
        //更新速度和位置
        _updatePosAndVel();
        
        //用于构建表面
        glm::vec3 tem=m_sphWallBox.min-glm::vec3(1.0);
        if(suf){
            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);
        }
    }

我们按照算法流程，首先求解领居结构，这在之前的章节里介绍过，这里就不多提了。

然后计算除了压力的所有外力，该函数_computerExternalForces()代码如下：

    void FluidSystem::_computerExternalForces(){
        float h2=m_smoothRadius*m_smoothRadius;//h^2
        for(int i=0;i<m_pointBuffer.size();i++){
            Point* pi=m_pointBuffer.get(i);
            pi->forces=glm::vec3(0.0);
            
            //邻居粒子装载
            pi->kernel_self=_computeNeighbor(i);
            m_neighborTable.point_commit();
            
            //外力计算
            int neighborCounts=m_neighborTable.getNeighborCounts(i);
            const float restVolume=m_pointMass/m_restDensity;
            for(unsigned 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;

                //F_Viscosity
                //m_kernelViscosity = 45.0f/(3.141592f * h^6);
                float vterm=m_kernelViscosity*m_viscosity*h_r*restVolume*restVolume;
                pi->forces-=(pi->velocity-pj->velocity)*vterm;
            }
            //F_gravity
            pi->forces+=m_gravityDir*m_pointMass;
            
            //F_boundary
            pi->forces+=_boundaryForce(pi)*m_pointMass;
            
            //初始化矫正因子
            pi->correction_pressure=0.0f;
            pi->correction_pressure_froce=glm::vec3(0.0);
        }
    }

计算完成外力后，进入矫正预测环节，该函数_predictionCorrection()代码为：

void FluidSystem::_predictionCorrection(){
    _density_error_too_large=true;
    int iteration=0;
    while ((iteration<minLoops)||((_density_error_too_large)&&(iteration<maxLoops))) {
        for(int i=0;i<m_pointBuffer.size();i++){
            Point* p=m_pointBuffer.get(i);
            _predictPositionAndVelocity(p);
        }
        
        //循环终止条件(所有粒子的最大预测密度)
        _max_predicted_density=0.0;
        
        for(int i=0;i<m_pointBuffer.size();i++){
            _computePredictedDensityAndPressure(i);
        }
        
        //循环终止条件
        float densityErrorInPercent=max(0.1f*_max_predicted_density-100.0f,0.0f);
        float maxDensityErrorAllowed=1.0f;
        
        //如果密度误差小于终止条件，则终止循环(小于密度误差阈值)
        if(densityErrorInPercent<maxDensityErrorAllowed)
            _density_error_too_large=false;
        
        for (int i=0; i<m_pointBuffer.size(); i++) {
            _computeCorrectivePressureForce(i);
        }
        iteration++;
    }
}

其中我们设置最小循环次数为3，最大循环次数为50。在矫正过程中，我们继续细分为3个步骤：

1).预测所有粒子的速度和位置，该函数_predictPositionAndVelocity(p)代码为：

void FluidSystem::_predictPositionAndVelocity(Point* p){
    // v' = v + delta_t * a
    // a = F / m
    // v' = v + delta_t * F * V / m
    // v' = v + delta_t * F * 1/density
    
    //计算预测速度和位置
    p->predicted_velocity=p->velocity+(p->forces+p->correction_pressure_froce)*m_deltaTime/m_pointMass;
    p->predicted_position=p->pos+p->predicted_velocity*m_deltaTime;
    
    //碰撞处理
    _collisionHandling(p->predicted_position,p->predicted_velocity);
    
    //初始化预测密度
    p->predicted_density=0.0;
}

值得注意的是，因为是预测的速度，所以需要好用矫正后的压力去计算。

2).计算预测的密度和压力，该函数_computePredictedDensityAndPressure(i)的代码为：

void FluidSystem::_computePredictedDensityAndPressure(int i){
    float h2=m_smoothRadius*m_smoothRadius;
    Point* pi=m_pointBuffer.get(i);
    int neighborCounts=m_neighborTable.getNeighborCounts(i);
    pi->predicted_density=pi->kernel_self*m_pointMass;
    //预测密度计算
    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);
        //需要用预测的位置去重新计算距离和核值
        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;
            pi->predicted_density+=m_kernelPoly6*pow(h2_r2, 3.0)*m_pointMass;
        }
    }
    auto sss=pi->predicted_density;
    //计算密度误差，仅考虑正向误差
    pi->density_error=max(pi->predicted_density-m_restDensity,0.0f);
    
    //更新压力,只矫正正压力，densityErrorFactor被预先计算并用作常数
    pi->correction_pressure+=max(pi->density_error*m_densityErrorFactor,0.0f);
    
    _max_predicted_density=max(_max_predicted_density,pi->predicted_density);
    
}

值得注意的是，我们这里都只考虑正向的误差，即压缩的密度误差和压强误差，拉伸的我们不做考虑。

3).计算矫正后的压力。该函数_computeCorrectivePressureForce(i)代码为：

void FluidSystem::_computeCorrectivePressureForce(int i){
    float h2=m_smoothRadius*m_smoothRadius;
    Point* pi=m_pointBuffer.get(i);
    int neighborCounts=m_neighborTable.getNeighborCounts(i);
    pi->correction_pressure_froce=glm::vec3(0.0f);
    
    float densSq=m_restDensity*m_restDensity;
    float pi_pres=pi->correction_pressure;
    
    //更新压力
    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);
        //需要用预测的位置去重新计算距离和核值
        glm::vec3 pi_pj=(pi->pos-pj->pos)*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){
            r=pow(r2,0.5);
            float h=m_smoothRadius;
            float pj_pres=pj->correction_pressure;
            glm::vec3 gradientKer=pi_pj*m_kernelSpiky/r*(h-r)*(h-r);
            float grad=pi_pres/densSq+pj_pres/(m_restDensity*m_restDensity);
            pi->correction_pressure_froce-=m_pointMass*m_pointMass*grad*gradientKer;
        }
    }
}

这里，我们的预测矫正迭代代码就结束了，我们推出迭代后，利用矫正好的压力求解新的速度和位置，该函数_updatePosAndVel()代码如下：

void FluidSystem::_updatePosAndVel(){
    for(int i=0;i<m_pointBuffer.size();i++){
        Point* pi=m_pointBuffer.get(i);
        pi->forces+=pi->correction_pressure_froce;
        const float invMass=1.0/m_pointMass;
        
//        //Symplectic Euler Integration
//        pi->velocity+=pi->forces*invMass*m_deltaTime;
//        pi->pos+=pi->velocity*m_deltaTime/m_unitScale;
        
        //Leapfrog integration
        glm::vec3 vnext = pi->velocity + pi->forces*invMass*m_deltaTime; // v(t+1/2) = v(t-1/2) + a(t) dt
        pi->velocity_eval = (pi->velocity + vnext)*0.5f;  //v(t+1) = [v(t-1/2) + v(t+1/2)] * 0.5
        pi->velocity = vnext;
        pi->pos += vnext*m_deltaTime/m_unitScale;  // 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;

    }
}

完成结果图：

对比普通SPH，我们发现在挤压处，PCISPH有明显的矫正。普通SPH效果如下：

0小龙虾0

关注

15
点赞
踩
31

收藏

觉得还不错? 一键收藏
15
评论
SPH（光滑粒子流体动力学）流体模拟实现五：PCISPH

SPH（光滑粒子流体动力学）流体模拟实现五：PCISPH我们知道真实的液体是不可压缩的，但我们在计算机中离散的计算流体运动，在一定的时间步长内，用标准的SPH方法求解，在在粒子聚集处容易发生挤压，造成压缩。有两种常用的方法模拟不可压缩性：1.在WCSPH(弱可压缩SPH)中，利用刚性状态方程（EOS）建模压力。2.通过求解压力泊松方程实现不可压缩性。但这两种方法都有很昂贵的计算费用。文章“Predictive-Corrective Incompressible SPH”中，提出了一种预测矫正的.
复制链接

扫一扫