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效果如下：

轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/203378.html

標籤：其他

上一篇：ubuntu 安裝 docker 并配置鏡像加速(使用 apt-get 進行安裝)

下一篇：Linux服務器安裝node