Poly6 kernel[1]

密度計算など.rは2乗項しか含まれていないので,平方根を計算しなくてよいのがメリット.

eq_poly6.gif

勾配は,

eq_poly6_gradient.gif

ラプラシアンは,

eq_poly6_laplacian.gif

グラフにすると,

kernel_poly6.png

ラプラシアンの値は第2軸(右側)参照

コード例

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
/*!
 * Poly6カーネル関数値の計算
 * @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カーネル関数勾配値の計算
 * @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カーネル関数ラプラシアンの計算
 * @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カーネルがr=0で勾配が0になり,パーティクルがクラスターを構成してしまう問題を解決する.

eq_spiky.gif

勾配は,

eq_spiky_gradient.gif

ラプラシアンは,

eq_spiky_laplacian.gif

[2]での定義は,

eq_spiky_desbrun.gif

グラフにすると,

kernel_spiky.png

ラプラシアンの値は第2軸(右側)参照

コード例

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
/*!
 * Spikyカーネル関数値の計算
 * @param[in] r 距離
 * @param[in] h 有効半径
 * @return 関数値
 */
inline double RxKernelSpiky(const double &r, const double &h)
{
    if(r >= 0.0 && r < h){
        double q = h-r;
        return 15/(RX_PI*pow(h, 6))*q*q*q;
    }
    else{
        return 0.0;
    }
}
 
/*!
 * Spikyカーネル関数勾配値の計算
 * @param[in] r 距離
 * @param[in] rij 相対位置ベクトル
 * @param[in] h 有効半径
 * @return 勾配値
 */
inline Vec3 RxKernelSpikyG(const double &r, const Vec3 &rij, const double &h)
{
    if(r >= 0.0 && r < h){
        double q = h-r;
        return -45.0/(RX_PI*pow(h, 6))*q*q*rij/r;
    }
    else{
        return Vec3(0.0);
    }
}
 
/*!
 * Spikyカーネル関数ラプラシアンの計算
 * @param[in] r 距離
 * @param[in] h 有効半径
 * @return ラプラシアンの値
 */
inline double RxKernelSpikyL(const double &r, const double &h)
{
    if(r >= 0.0 && r < h){
        double q = h-r;
        return -90.0/(RX_PI*pow(h, 6))*(q*q/r-q);
    }
    else{
        return 0.0;
    }
}

hが定数ならば,上式のα部分はグローバル変数にでもあらかじめ計算して入れておいて使った方がよい.

Viscosity kernel[1]

粘性拡散項計算用.Laplacianが線型な関数になる.

eq_visc.gif

勾配は,

eq_visc_gradient.gif

ラプラシアンは,

eq_visc_laplacian.gif

グラフにすると,

kernel_viscosity.png

ラプラシアンの値は第2軸(右側)参照

コード例

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
/*!
 * Viscosityカーネル関数値の計算
 * @param[in] r 距離
 * @param[in] h 有効半径
 * @return 関数値
 */
inline double RxKernelViscosity(const double &r, const double &h)
{
    if(r >= 0.0 && r < h){
        double q = r/h;
        return 15.0/(2.0*RX_PI*pow(h, 3))*(-0.5*q*q*q+q*q+0.5/q-1);
    }
    else{
        return 0.0;
    }
}
 
/*!
 * Viscosityカーネル関数勾配値の計算
 * @param[in] r 距離
 * @param[in] rij 相対位置ベクトル
 * @param[in] h 有効半径
 * @return 勾配値
 */
inline Vec3 RxKernelViscosityG(const double &r, const Vec3 &rij, const double &h)
{
    if(r >= 0.0 && r < h){
        double q = r/h;
        return 15.0/(2.0*RX_PI*pow(h, 5))*(-1.5*q+2-0.5/(q*q*q))*rij;
    }
    else{
        return Vec3(0.0);
    }
}
 
/*!
 * Viscosityカーネル関数ラプラシアンの計算
 * @param[in] r 距離
 * @param[in] h 有効半径
 * @return ラプラシアンの値
 */
inline double RxKernelViscosityL(const double &r, const double &h)
{
    if(r >= 0.0 && r < h){
        double q = h-r;
        return 45.0/(RX_PI*pow(h, 6))*q;
    }
    else{
        return 0.0;
    }
}

hが定数ならば,上式のα部分はグローバル変数にでもあらかじめ計算して入れておいて使った方がよい.

Spline kernel[3,4]

sph_kernel_spline.eq1.gif

ここで,

sph_kernel_spline.eq2.gif

勾配は,

sph_kernel_spline.eq3.gif

ここで,

sph_kernel_spline.eq4.gif

ラプラシアンは,

sph_kernel_spline.eq5.gif

ここで,

sph_kernel_spline.eq6.gif

グラフにすると,

kernel_spline.png

ラプラシアンの値は第2軸(右側)参照

Super Gaussian kernel[3]

eq_gaussian2.gif

勾配は,

eq_gaussian2_gradient.gif

ラプラシアンは,

eq_gaussian2_laplacian.gif

グラフにすると,

kernel_gaussian2.png

ラプラシアンの値は第2軸(右側)参照

コード例

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
/*!
 * Gaussianカーネル関数値の計算
 * @param[in] r 距離
 * @param[in] h 有効半径
 * @return 関数値
 */
inline double RxKernelSuperGaussian(const double &r, const double &h)
{
    double q = r/h;
    return 1.0/(pow(RX_PI, 1.5)*pow(h, 3))*(2.5-r*r)*exp(-q*q);
    //return 1.0/(sqrt(RX_PI)*h)*exp(-q*q);
}
 
/*!
 * Gaussianカーネル関数勾配値の計算
 * @param[in] r 距離
 * @param[in] rij 相対位置ベクトル
 * @param[in] h 有効半径
 * @return 勾配値
 */
inline Vec3 RxKernelSuperGaussianG(const double &r, const Vec3 &rij, const double &h)
{
    double q = r/h;
    return -2.0/(pow(RX_PI, 1.5)*pow(h, 3))*(1.0+2.5/(h*h)-q*q)*exp(-q*q)*rij;
}
 
/*!
 * Gaussianカーネル関数ラプラシアンの計算
 * @param[in] r 距離
 * @param[in] h 有効半径
 * @return ラプラシアンの値
 */
inline double RxKernelSuperGaussianL(const double &r, const double &h)
{
    double q = r/h;
    return 2.0/(pow(RX_PI, 1.5)*pow(h, 3))*(2.0*q*q-(1.0+2.5/(h*h)-q*q)*(1.0-2.0*q*q))*exp(-q*q);
}

hが定数ならば,上式のα部分はグローバル変数にでもあらかじめ計算して入れておいて使った方がよい.

Fourth-order weighting function

eq_fourth.gif

勾配は,

eq_fourth_gradient.gif

ラプラシアンは,

eq_fourth_laplacian.gif

グラフにすると,

kernel_fourth.png

ラプラシアンの値は第2軸(右側)参照

コード例

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
/*!
 * Fourth-order重み関数値の計算
 * @param[in] r 距離
 * @param[in] h 有効半径
 * @return 関数値
 */
inline double RxKernelFourth(const double &r, const double &h)
{
    double q = 3*r/h;
    if(q >= 0.0 && q < 1.0){
        double q3 = (3-q)*(3-q)*(3-q)*(3-q)*(3-q);
        double q2 = (2-q)*(2-q)*(2-q)*(2-q)*(2-q);
        double q1 = (1-q)*(1-q)*(1-q)*(1-q)*(1-q);
        return 9.0/(40.0*RX_PI*pow(h, 3))*(q3-6.0*q2+15.0*q1);
    }
    else if(q >= 1.0 && q < 2.0){
        double q3 = (3-q)*(3-q)*(3-q)*(3-q)*(3-q);
        double q2 = (2-q)*(2-q)*(2-q)*(2-q)*(2-q);
        return 9.0/(40.0*RX_PI*pow(h, 3))*(q3-6.0*q2);
    }
    else if(q >= 2.0 && q < 3.0){
        double q3 = (3-q)*(3-q)*(3-q)*(3-q)*(3-q);
        return 9.0/(40.0*RX_PI*pow(h, 3))*(q3);
    }
    else{
        return 0.0;
    }
}
 
/*!
 * Fourth-order重み関数勾配値の計算
 * @param[in] r 距離
 * @param[in] rij 相対位置ベクトル
 * @param[in] h 有効半径
 * @return 勾配値
 */
inline Vec3 RxKernelFourthG(const double &r, const Vec3 &rij, const double &h)
{
    double q = 3*r/h;
    if(q >= 0.0 && q < 1.0){
        double q3 = (3-q)*(3-q)*(3-q)*(3-q);
        double q2 = (2-q)*(2-q)*(2-q)*(2-q);
        double q1 = (1-q)*(1-q)*(1-q)*(1-q);
        return -27.0/(8.0*RX_PI*pow(h, 4))*(q3-6.0*q2+15.0*q1)*rij/r;
    }
    else if(q >= 1.0 && q < 2.0){
        double q3 = (3-q)*(3-q)*(3-q)*(3-q);
        double q2 = (2-q)*(2-q)*(2-q)*(2-q);
        return -27.0/(8.0*RX_PI*pow(h, 4))*(q3-6.0*q2)*rij/r;
    }
    else if(q >= 2.0 && q < 3.0){
        double q3 = (3-q)*(3-q)*(3-q)*(3-q);
        return -27.0/(8.0*RX_PI*pow(h, 4))*(q3)*rij/r;
    }
    else{
        return Vec3(0.0);
    }
}
 
/*!
 * Fourth-order重み関数ラプラシアンの計算
 * @param[in] r 距離
 * @param[in] h 有効半径
 * @return ラプラシアンの値
 */
inline double RxKernelFourthL(const double &r, const double &h)
{
    double q = 3*r/h;
    if(q >= 0.0 && q < 1.0){
        double q3 = (3-q)*(3-q)*(3-q);
        double q2 = (2-q)*(2-q)*(2-q);
        double q1 = (1-q)*(1-q)*(1-q);
        return 81.0/(2.0*RX_PI*pow(h, 5))*(q3-6.0*q2+15.0*q1);
    }
    else if(q >= 1.0 && q < 2.0){
        double q3 = (3-q)*(3-q)*(3-q);
        double q2 = (2-q)*(2-q)*(2-q);
        return 81.0/(2.0*RX_PI*pow(h, 5))*(q3-6.0*q2);
    }
    else if(q >= 2.0 && q < 3.0){
        double q3 = (3-q)*(3-q)*(3-q);
        return 81.0/(2.0*RX_PI*pow(h, 5))*(q3);
    }
    else{
        return 0.0;
    }
}

hが定数ならば,上式のα部分はグローバル変数にでもあらかじめ計算して入れておいて使った方がよい.

カーネル関数の係数計算

SPH法ではカーネル関数Wは以下の条件を満たすことを仮定している.

sph_kernel.eq1.gif

これを満たすように各カーネルの係数が決定される. 例としてPoly6カーネル(sph_kernel.eq2.gif)の場合で実際にやってみよう.

sph_kernel.eq3.gif

ここで係数をsph_kernel.eq4.gifとした.

2次元では有効半径hの円を使って積分する. 下の図のように微少範囲sph_kernel.eq5.gifを積分する.

sph_kernel.eq13.gif

これが1となるためには

sph_kernel.eq6.gif

と導かれる.

kernel2d.jpg
2次元の場合

次に3次元の場合を考えてみよう. 下図左のように球の中心から表面まで伸びた角錐形状内の微少体積を考える. 各辺の長さはsph_kernel.eq7.gif, sph_kernel.eq8.gif, sph_kernel.eq9.gif(図右参照)となるので, カーネル関数の積分は以下のように表される.

sph_kernel.eq10.gif

Poly6カーネルを代入して計算する.

sph_kernel.eq14.gif

sph_kernel.eq11.gifとなるためには

sph_kernel.eq12.gif

と導かれる.

kernel3d.jpg
3次元の場合

参考文献

[1] Muller et al., Particle-Based Fluid Simulation for Interactive Applications, SCA2003.

[2] Desbrun and Cani, Smoothed Particles: A new paradigm for animating highly deformable bodies, Eurographics Workshop on Computer Animation and Simulation (EGCAS), 1996.

[3] Monaghan, Smoothed particle hydrodynamics, Annual Review of Astronomy and Astrophysics, 30, pp543-574, 1992.

[4] Becker and Teschner, Weakly compressible SPH for free surface flows, SCA2007.


添付ファイル: filesph_kernel.eq14.gif 620件 [詳細] filesph_kernel.eq13.gif 597件 [詳細] filekernel3d.jpg 1346件 [詳細] filekernel2d.jpg 1483件 [詳細] filesph_kernel.eq10.gif 549件 [詳細] filesph_kernel.eq3.gif 505件 [詳細] filesph_kernel.eq12.gif 609件 [詳細] filesph_kernel.eq1.gif 617件 [詳細] filesph_kernel.eq9.gif 553件 [詳細] filesph_kernel.eq4.gif 551件 [詳細] filesph_kernel.eq2.gif 549件 [詳細] filesph_kernel.eq8.gif 552件 [詳細] filesph_kernel.eq5.gif 551件 [詳細] filesph_kernel.eq6.gif 589件 [詳細] filesph_kernel.eq11.gif 534件 [詳細] filesph_kernel.eq7.gif 592件 [詳細] filesph_kernel_spline.eq5.gif 656件 [詳細] filesph_kernel_spline.eq2.gif 620件 [詳細] filesph_kernel_spline.eq6.gif 664件 [詳細] filesph_kernel_spline.eq3.gif 635件 [詳細] filesph_kernel_spline.eq4.gif 639件 [詳細] filesph_kernel_spline.eq1.gif 644件 [詳細] filekernel_spline.png 1822件 [詳細] filekernel_gaussian2.png 2185件 [詳細] filekernel_viscosity.png 1624件 [詳細] filekernel_poly6.png 1866件 [詳細] fileeq_visc_laplacian.gif 1745件 [詳細] filekernel_spiky.png 1536件 [詳細] filekernel_fourth.png 1663件 [詳細] fileeq_visc.gif 1668件 [詳細] fileeq_visc_gradient.gif 1683件 [詳細] fileeq_spline_gradient.gif 737件 [詳細] fileeq_spiky_gradient.gif 1578件 [詳細] fileeq_spiky_desbrun.gif 1684件 [詳細] fileeq_spiky_laplacian.gif 1693件 [詳細] fileeq_gaussian2.gif 1576件 [詳細] fileeq_poly6.gif 2066件 [詳細] fileeq_gaussian2_gradient.gif 1548件 [詳細] fileeq_gaussian2_laplacian.gif 1647件 [詳細] fileeq_poly6_gradient.gif 1808件 [詳細] fileeq_poly6_laplacian.gif 1611件 [詳細] fileeq_spiky.gif 1706件 [詳細] fileeq_gaussian1.gif 733件 [詳細] fileeq_fourth.gif 1740件 [詳細] fileeq_fourth_laplacian.gif 1719件 [詳細] fileeq_fourth_gradient.gif 1680件 [詳細]

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2014-02-06 (木) 17:37:12 (2418d)