SPH法の重み関数
をテンプレートにして作成
[
トップ
|
新規
|
一覧
|
検索
|
最終更新
|
ヘルプ
]
開始行:
----
#contents
----
***Poly6 kernel[1] [#hcde5047]
密度計算など．rは2乗項しか含まれていないので，平方根を計...
#ref(eq_poly6.gif)
勾配は，
#ref(eq_poly6_gradient.gif)
ラプラシアンは，
#ref(eq_poly6_laplacian.gif)
グラフにすると，
#ref(kernel_poly6.png,,50%)
ラプラシアンの値は第2軸(右側)参照
コード例
#code(C){{
/*!
* Poly6カーネル関数値の計算
* @param[in] r 距離
* @param[in] h 有効半径
* @return 関数値
*/
inline double RxKernelPoly6(const double &r, const double...
{
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 &r...
{
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 doubl...
{
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] [#yb8b47cb]
圧力項計算用．Poly6カーネルがr=0で勾配が0になり，パーティ...
#ref(eq_spiky.gif)
勾配は，
#ref(eq_spiky_gradient.gif)
ラプラシアンは，
#ref(eq_spiky_laplacian.gif)
[2]での定義は，
#ref(eq_spiky_desbrun.gif)
グラフにすると，
#ref(kernel_spiky.png,,50%)
ラプラシアンの値は第2軸(右側)参照
コード例
#code(C){{
/*!
* Spikyカーネル関数値の計算
* @param[in] r 距離
* @param[in] h 有効半径
* @return 関数値
*/
inline double RxKernelSpiky(const double &r, const double...
{
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 &r...
{
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 doubl...
{
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] [#a85e448a]
粘性拡散項計算用．Laplacianが線型な関数になる．
#ref(eq_visc.gif)
勾配は，
#ref(eq_visc_gradient.gif)
ラプラシアンは，
#ref(eq_visc_laplacian.gif)
グラフにすると，
#ref(kernel_viscosity.png,,50%)
ラプラシアンの値は第2軸(右側)参照
コード例
#code(C){{
/*!
* Viscosityカーネル関数値の計算
* @param[in] r 距離
* @param[in] h 有効半径
* @return 関数値
*/
inline double RxKernelViscosity(const double &r, const do...
{
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...
}
else{
return 0.0;
}
}
/*!
* Viscosityカーネル関数勾配値の計算
* @param[in] r 距離
* @param[in] rij 相対位置ベクトル
* @param[in] h 有効半径
* @return 勾配値
*/
inline Vec3 RxKernelViscosityG(const double &r, const Vec...
{
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)...
}
else{
return Vec3(0.0);
}
}
/*!
* Viscosityカーネル関数ラプラシアンの計算
* @param[in] r 距離
* @param[in] h 有効半径
* @return ラプラシアンの値
*/
inline double RxKernelViscosityL(const double &r, const d...
{
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] [#o437ed20]
#ref(sph_kernel_spline.eq1.gif,nolink,70%)
ここで，
#ref(sph_kernel_spline.eq2.gif,nolink,70%)
勾配は，
#ref(sph_kernel_spline.eq3.gif,nolink,70%)
ここで，
#ref(sph_kernel_spline.eq4.gif,nolink,70%)
ラプラシアンは，
#ref(sph_kernel_spline.eq5.gif,nolink,70%)
ここで，
#ref(sph_kernel_spline.eq6.gif,nolink,70%)
グラフにすると，
#ref(kernel_spline.png,,50%)
ラプラシアンの値は第2軸(右側)参照
***Super Gaussian kernel[3] [#s1074f09]
#ref(eq_gaussian2.gif)
勾配は，
#ref(eq_gaussian2_gradient.gif)
ラプラシアンは，
#ref(eq_gaussian2_laplacian.gif)
グラフにすると，
#ref(kernel_gaussian2.png,,50%)
ラプラシアンの値は第2軸(右側)参照
コード例
#code(C){{
/*!
* Gaussianカーネル関数値の計算
* @param[in] r 距離
* @param[in] h 有効半径
* @return 関数値
*/
inline double RxKernelSuperGaussian(const double &r, cons...
{
double q = r/h;
return 1.0/(pow(RX_PI, 1.5)*pow(h, 3))*(2.5-r*r)*exp(-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...
{
double q = r/h;
return -2.0/(pow(RX_PI, 1.5)*pow(h, 3))*(1.0+2.5/(h*h)-q...
}
/*!
* Gaussianカーネル関数ラプラシアンの計算
* @param[in] r 距離
* @param[in] h 有効半径
* @return ラプラシアンの値
*/
inline double RxKernelSuperGaussianL(const double &r, con...
{
double q = r/h;
return 2.0/(pow(RX_PI, 1.5)*pow(h, 3))*(2.0*q*q-(1.0+2.5...
}
}}
hが定数ならば，上式のα部分はグローバル変数にでもあらかじ...
***Fourth-order weighting function [#u853dbcc]
#ref(eq_fourth.gif)
勾配は，
#ref(eq_fourth_gradient.gif)
ラプラシアンは，
#ref(eq_fourth_laplacian.gif)
グラフにすると，
#ref(kernel_fourth.png,,50%)
ラプラシアンの値は第2軸(右側)参照
コード例
#code(C){{
/*!
* Fourth-order重み関数値の計算
* @param[in] r 距離
* @param[in] h 有効半径
* @return 関数値
*/
inline double RxKernelFourth(const double &r, const doubl...
{
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 &...
{
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)*...
}
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 doub...
{
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が定数ならば，上式のα部分はグローバル変数にでもあらかじ...
**カーネル関数の係数計算 [#y6434ffe]
SPH法ではカーネル関数Wは以下の条件を満たすことを仮定して...
#ref(sph_kernel.eq1.gif,nolink,70%)
これを満たすように各カーネルの係数が決定される．
例としてPoly6カーネル(&ref(sph_kernel.eq2.gif,nolink,70%)...
#ref(sph_kernel.eq3.gif,nolink,70%)
ここで係数を&ref(sph_kernel.eq4.gif,nolink,70%);とした．
2次元では有効半径hの円を使って積分する．
下の図のように微少範囲&ref(sph_kernel.eq5.gif,nolink,70%)...
#ref(sph_kernel.eq13.gif,nolink,70%)
これが1となるためには
#ref(sph_kernel.eq6.gif,nolink,70%)
と導かれる．
|&ref(kernel2d.jpg,,100%);|
|CENTER:2次元の場合|
次に3次元の場合を考えてみよう．
下図左のように球の中心から表面まで伸びた角錐形状内の微少...
各辺の長さは&ref(sph_kernel.eq7.gif,nolink,70%);, &ref(sp...
カーネル関数の積分は以下のように表される．
#ref(sph_kernel.eq10.gif,nolink,70%)
Poly6カーネルを代入して計算する．
#ref(sph_kernel.eq14.gif,nolink,70%)
&ref(sph_kernel.eq11.gif,nolink,70%);となるためには
#ref(sph_kernel.eq12.gif,nolink,70%)
と導かれる．
|&ref(kernel3d.jpg,,100%);|
|CENTER:3次元の場合|
**参考文献 [#zeacb0e2]
>> [1] Muller et al., Particle-Based Fluid Simulation for...
>> [2] Desbrun and Cani, Smoothed Particles: A new paradi...
>> [3] Monaghan, Smoothed particle hydrodynamics, Annual ...
>> [4] Becker and Teschner, Weakly compressible SPH for f...
終了行:
----
#contents
----
***Poly6 kernel[1] [#hcde5047]
密度計算など．rは2乗項しか含まれていないので，平方根を計...
#ref(eq_poly6.gif)
勾配は，
#ref(eq_poly6_gradient.gif)
ラプラシアンは，
#ref(eq_poly6_laplacian.gif)
グラフにすると，
#ref(kernel_poly6.png,,50%)
ラプラシアンの値は第2軸(右側)参照
コード例
#code(C){{
/*!
* Poly6カーネル関数値の計算
* @param[in] r 距離
* @param[in] h 有効半径
* @return 関数値
*/
inline double RxKernelPoly6(const double &r, const double...
{
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 &r...
{
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 doubl...
{
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] [#yb8b47cb]
圧力項計算用．Poly6カーネルがr=0で勾配が0になり，パーティ...
#ref(eq_spiky.gif)
勾配は，
#ref(eq_spiky_gradient.gif)
ラプラシアンは，
#ref(eq_spiky_laplacian.gif)
[2]での定義は，
#ref(eq_spiky_desbrun.gif)
グラフにすると，
#ref(kernel_spiky.png,,50%)
ラプラシアンの値は第2軸(右側)参照
コード例
#code(C){{
/*!
* Spikyカーネル関数値の計算
* @param[in] r 距離
* @param[in] h 有効半径
* @return 関数値
*/
inline double RxKernelSpiky(const double &r, const double...
{
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 &r...
{
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 doubl...
{
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] [#a85e448a]
粘性拡散項計算用．Laplacianが線型な関数になる．
#ref(eq_visc.gif)
勾配は，
#ref(eq_visc_gradient.gif)
ラプラシアンは，
#ref(eq_visc_laplacian.gif)
グラフにすると，
#ref(kernel_viscosity.png,,50%)
ラプラシアンの値は第2軸(右側)参照
コード例
#code(C){{
/*!
* Viscosityカーネル関数値の計算
* @param[in] r 距離
* @param[in] h 有効半径
* @return 関数値
*/
inline double RxKernelViscosity(const double &r, const do...
{
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...
}
else{
return 0.0;
}
}
/*!
* Viscosityカーネル関数勾配値の計算
* @param[in] r 距離
* @param[in] rij 相対位置ベクトル
* @param[in] h 有効半径
* @return 勾配値
*/
inline Vec3 RxKernelViscosityG(const double &r, const Vec...
{
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)...
}
else{
return Vec3(0.0);
}
}
/*!
* Viscosityカーネル関数ラプラシアンの計算
* @param[in] r 距離
* @param[in] h 有効半径
* @return ラプラシアンの値
*/
inline double RxKernelViscosityL(const double &r, const d...
{
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] [#o437ed20]
#ref(sph_kernel_spline.eq1.gif,nolink,70%)
ここで，
#ref(sph_kernel_spline.eq2.gif,nolink,70%)
勾配は，
#ref(sph_kernel_spline.eq3.gif,nolink,70%)
ここで，
#ref(sph_kernel_spline.eq4.gif,nolink,70%)
ラプラシアンは，
#ref(sph_kernel_spline.eq5.gif,nolink,70%)
ここで，
#ref(sph_kernel_spline.eq6.gif,nolink,70%)
グラフにすると，
#ref(kernel_spline.png,,50%)
ラプラシアンの値は第2軸(右側)参照
***Super Gaussian kernel[3] [#s1074f09]
#ref(eq_gaussian2.gif)
勾配は，
#ref(eq_gaussian2_gradient.gif)
ラプラシアンは，
#ref(eq_gaussian2_laplacian.gif)
グラフにすると，
#ref(kernel_gaussian2.png,,50%)
ラプラシアンの値は第2軸(右側)参照
コード例
#code(C){{
/*!
* Gaussianカーネル関数値の計算
* @param[in] r 距離
* @param[in] h 有効半径
* @return 関数値
*/
inline double RxKernelSuperGaussian(const double &r, cons...
{
double q = r/h;
return 1.0/(pow(RX_PI, 1.5)*pow(h, 3))*(2.5-r*r)*exp(-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...
{
double q = r/h;
return -2.0/(pow(RX_PI, 1.5)*pow(h, 3))*(1.0+2.5/(h*h)-q...
}
/*!
* Gaussianカーネル関数ラプラシアンの計算
* @param[in] r 距離
* @param[in] h 有効半径
* @return ラプラシアンの値
*/
inline double RxKernelSuperGaussianL(const double &r, con...
{
double q = r/h;
return 2.0/(pow(RX_PI, 1.5)*pow(h, 3))*(2.0*q*q-(1.0+2.5...
}
}}
hが定数ならば，上式のα部分はグローバル変数にでもあらかじ...
***Fourth-order weighting function [#u853dbcc]
#ref(eq_fourth.gif)
勾配は，
#ref(eq_fourth_gradient.gif)
ラプラシアンは，
#ref(eq_fourth_laplacian.gif)
グラフにすると，
#ref(kernel_fourth.png,,50%)
ラプラシアンの値は第2軸(右側)参照
コード例
#code(C){{
/*!
* Fourth-order重み関数値の計算
* @param[in] r 距離
* @param[in] h 有効半径
* @return 関数値
*/
inline double RxKernelFourth(const double &r, const doubl...
{
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 &...
{
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)*...
}
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 doub...
{
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が定数ならば，上式のα部分はグローバル変数にでもあらかじ...
**カーネル関数の係数計算 [#y6434ffe]
SPH法ではカーネル関数Wは以下の条件を満たすことを仮定して...
#ref(sph_kernel.eq1.gif,nolink,70%)
これを満たすように各カーネルの係数が決定される．
例としてPoly6カーネル(&ref(sph_kernel.eq2.gif,nolink,70%)...
#ref(sph_kernel.eq3.gif,nolink,70%)
ここで係数を&ref(sph_kernel.eq4.gif,nolink,70%);とした．
2次元では有効半径hの円を使って積分する．
下の図のように微少範囲&ref(sph_kernel.eq5.gif,nolink,70%)...
#ref(sph_kernel.eq13.gif,nolink,70%)
これが1となるためには
#ref(sph_kernel.eq6.gif,nolink,70%)
と導かれる．
|&ref(kernel2d.jpg,,100%);|
|CENTER:2次元の場合|
次に3次元の場合を考えてみよう．
下図左のように球の中心から表面まで伸びた角錐形状内の微少...
各辺の長さは&ref(sph_kernel.eq7.gif,nolink,70%);, &ref(sp...
カーネル関数の積分は以下のように表される．
#ref(sph_kernel.eq10.gif,nolink,70%)
Poly6カーネルを代入して計算する．
#ref(sph_kernel.eq14.gif,nolink,70%)
&ref(sph_kernel.eq11.gif,nolink,70%);となるためには
#ref(sph_kernel.eq12.gif,nolink,70%)
と導かれる．
|&ref(kernel3d.jpg,,100%);|
|CENTER:3次元の場合|
**参考文献 [#zeacb0e2]
>> [1] Muller et al., Particle-Based Fluid Simulation for...
>> [2] Desbrun and Cani, Smoothed Particles: A new paradi...
>> [3] Monaghan, Smoothed particle hydrodynamics, Annual ...
>> [4] Becker and Teschner, Weakly compressible SPH for f...
ページ名:
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
