【CG物理模拟系列】流体模拟--粒子法之SPH法的加权函数计算

SPH法中的加权函数

Poly6 kernel[1]

用于密度计算等。由于r只存在于2次项中，可以省去平方根的计算。

梯度为，

拉普拉斯算子为，

如图所示，

右侧坐标为拉普拉斯值

Code Sample

  /*!

 * Poly6 kernel函数值的计算
 * @param[in] r 距离
 * @param[in] h 有效半径
 * @return 函数值
 */inline double RxKernelPoly6(const double &r, const double &h)
{
    if(r >= 0.0 && r < h){
        double q = h*h-r*r;
        return 315.0/(64.0*RX_PI*pow(h, 9))*q*q*q;
    }
    else{
        return 0.0;
    }
}
 
/*!
 * Poly6 kernel梯度值的计算
 * @param[in] r 距离
 * @param[in] rij 相对位置
 * @param[in] h 有效半径
 * @return 梯度值
 */inline Vec3 RxKernelPoly6G(const double &r, const Vec3 &rij, const double &h)
{
    if(r >= 0.0 && r < h){
        double q = h*h-r*r;
        return -945.0/(32.0*RX_PI*pow(h, 9))*q*q*rij;
    }
    else{
        return Vec3(0.0);
    }
}
 
/*!
 * Poly6 kernel函数拉普拉斯算子的计算
 * @param[in] r 距离
 * @param[in] h 有效半径
 * @return 拉普拉斯值

 */ inline  double RxKernelPoly6L( const  double &r,  const  double &h)

{
    if(r >= 0.0 && r < h){
        double q = h*h-r*r;
        return -945.0/(32.0*RX_PI*pow(h, 9))*(3*q*q-4*r*r*q);
    }
    else{
        return 0.0;
    }
}

h为定值的时候，上式的α部分作为全局变量计算即可。

↑

Spiky kernel[1,2] 

用于计算圧力项。当Poly6 kernel的r=0时，其梯度也为0，把粒子当作是由cluster构成的问题求解即可。

梯度为，

拉普拉斯值为，

[2]的定义为，