(57)【特許請求の範囲】
【請求項１】
流体の経時的な形状変化をシミュレーションする場合に、不自然な分布の速度分布が発生する挙動を抑制するシミュレーションプログラムにおいて、
コンピュータに、
格納部から、前記流体を形成する粒子の位置ベクトル及び速度ベクトルを取得し、
前記流体を粒子の集合で表現した流体モデルについて、取得した前記位置ベクトル及び速度ベクトルを用いて、前記流体を形成する各粒子に係る力および該各粒子に係る表面張力を算出し、
前記速度ベクトルに算出した各粒子に係る力及び各粒子に係る表面張力を加えて該各粒子の速度を更新し、
更新した該各粒子の速度に基づいて、前記流体の形状変化を算出し、
算出された前記形状変化を出力部に出力する
処理を実行させるシミュレーションプログラムであって、
前記該各粒子の速度の更新においては、
【数６５】
0005842992.tif 000147
を用い、前記各粒子に係る力は、上記（式１２）の右辺第２項で示される圧力勾配項、右辺第３項で示される粘性応力項、右辺第４項で示される重力項、右辺第５項は粒子間距離が近づきすぎないようにするための力であり、右辺第６項は表面張力であり、ここで、
【数６６】
0005842992.tif 000148
【数６７】
0005842992.tif 000149
【数６８】
0005842992.tif 000150
また、下付き添え字は粒子の番号を表し、
0005842992.tif 000151
：i番目の粒子の位置ベクトル、速度ベクトル、密度、圧力であり、
0005842992.tif 000152
：j番目の質量であり、
0005842992.tif 000153
：粘性項を計算するために導入した定数であり、
【数６９】
0005842992.tif 000154
であり
0005842992.tif 000155
【数７０】
0005842992.tif 000156
ｙ＝０の平面が固体の壁面の場合、
0005842992.tif 000157
であり、
ｇ：重力加速度
ｙ：ｙ方向の単位ベクトル
ｈ：カーネル関数の値が非ゼロとなる球形領域の半径である影響半径
Ｗ：ＳＰＨ法で用いられるカーネル関数の重ね合わせ
0005842992.tif 000158
は
【数７１】
0005842992.tif 000159
の再標準化行列であり、
【数７２】
0005842992.tif 000160
とおくと、
0005842992.tif 000161
0005842992.tif 000162
であり、
0005842992.tif 000163
は粒子ｉ，ｊ間での１次元リーマン問題を解いて求めた時空間の中間値であり、以下のように決定し、
まず、粒子ｉ，ｊに関して、以下の特性量を定め、
【数７３】
0005842992.tif 000164
【数７４】
0005842992.tif 000165
【数７５】
0005842992.tif 000166
【数７６】
0005842992.tif 000167
さらに、各勾配を以下のよう計算し、
【数７７】
0005842992.tif 000168
【数７８】
0005842992.tif 000169
【数７９】
0005842992.tif 000170
【数８０】
0005842992.tif 000171
【数８１】
0005842992.tif 000172
【数８２】
0005842992.tif 000173
ここで、
0005842992.tif 000174
0005842992.tif 000175
とし、上記の量を用いて、
0005842992.tif 000176
は以下のように決定する
【数８３】
0005842992.tif 000177
【数８４】
0005842992.tif 000178
【数８５】
0005842992.tif 000179
【数８６】
0005842992.tif 000180
【数８７】
0005842992.tif 000181
ことを特徴とするシミュレーションプログラム。
【請求項２】
前記表面張力の算出は、空気との界面の表面エネルギーを算出した結果と前記空気以外との界面の表面エネルギーを算出した結果の和を算出することを特徴とする請求項１記載のシミュレーションプログラム。
【請求項３】
前記形状変化の算出は、前記流体が他相を通り抜けようとする場合、前記表面張力の算出において用いられた前記粒子の位置及び速度を補正することを特徴とする請求項１記載のシミュレーションプログラム。
【請求項４】
流体の経時的な形状変化をシミュレーションする場合に、不自然な分布の速度分布が発生する挙動を抑制するシミュレーション方法において、
コンピュータが、
格納部から、前記流体を形成する粒子の位置ベクトル及び速度ベクトルを取得し、
前記流体を粒子の集合で表現した流体モデルについて、取得した前記位置ベクトル及び速度ベクトルを用いて、前記流体を形成する各粒子に係る力および該各粒子に係る表面張力を算出し、
前記速度ベクトルに算出した各粒子に係る力及び各粒子に係る表面張力を加えて該各粒子の速度を更新し、
更新した該各粒子の速度に基づいて、前記流体の形状変化を算出し、
算出された前記形状変化を出力部に出力する
シミュレーション方法であって、
前記該各粒子の速度の更新においては、
【数８８】
0005842992.tif 000182
を用い、前記各粒子に係る力は、上記（式１２）の右辺第２項で示される圧力勾配項、右辺第３項で示される粘性応力項、右辺第４項で示される重力項、右辺第５項は粒子間距離が近づきすぎないようにするための力であり、右辺第６項は表面張力であり、ここで、
【数８９】
0005842992.tif 000183
【数９０】
0005842992.tif 000184
【数９１】
0005842992.tif 000185
また、下付き添え字は粒子の番号を表し、
0005842992.tif 000186
：i番目の粒子の位置ベクトル、速度ベクトル、密度、圧力であり、
0005842992.tif 000187
：j番目の質量であり、
0005842992.tif 000188
：粘性項を計算するために導入した定数であり、
【数９２】
0005842992.tif 000189
であり
0005842992.tif 000190
【数９３】
0005842992.tif 000191
ｙ＝０の平面が固体の壁面の場合、
0005842992.tif 000192
であり、
ｇ：重力加速度
ｙ：ｙ方向の単位ベクトル
ｈ：カーネル関数の値が非ゼロとなる球形領域の半径である影響半径
Ｗ：ＳＰＨ法で用いられるカーネル関数の重ね合わせ
0005842992.tif 000193
は
【数９４】
0005842992.tif 000194
の再標準化行列であり、
【数９５】
0005842992.tif 000195
とおくと、
0005842992.tif 000196
0005842992.tif 000197
であり、
0005842992.tif 000198
は粒子ｉ，ｊ間での１次元リーマン問題を解いて求めた時空間の中間値であり、以下のように決定し、
まず、粒子ｉ，ｊに関して、以下の特性量を定め、
【数９６】
0005842992.tif 000199
【数９７】
0005842992.tif 000200
【数９８】
0005842992.tif 000201
【数９９】
0005842992.tif 000202
さらに、各勾配を以下のよう計算し、
【数１００】
0005842992.tif 000203
【数１０１】
0005842992.tif 000204
【数１０２】
0005842992.tif 000205
【数１０３】
0005842992.tif 000206
【数１０４】
0005842992.tif 000207
【数１０５】
0005842992.tif 000208
ここで、
0005842992.tif 000209
0005842992.tif 000210
とし、上記の量を用いて、
0005842992.tif 000211
は以下のように決定する
【数１０６】
0005842992.tif 000212
【数１０７】
0005842992.tif 000213
【数１０８】
0005842992.tif 000214
【数１０９】
0005842992.tif 000215
【数１１０】
0005842992.tif 000216
ことを特徴とするシミュレーション方法。
【請求項５】
流体の経時的な形状変化をシミュレーションする場合に、不自然な分布の速度分布が発生する挙動を抑制するシミュレーション装置において、
前記流体を形成する粒子の位置ベクトル及び速度ベクトルを格納する格納部と、
前記流体を粒子の集合で表現した流体モデルについて、前記位置ベクトル及び速度ベクトルを用いて、前記流体を形成する各粒子に係る力および該各粒子に係る表面張力を算出
する表面張力算出部と、
前記速度ベクトルに算出した各粒子に係る力及び各粒子に係る表面張力を加えて該各粒子の速度を更新する更新部と、
更新した該各粒子の速度に基づいて、前記流体の形状変化を算出する形状変化算出部と、
算出された前記形状変化を出力する出力部と、
を備えるシミュレーション装置であって、
前記該各粒子の速度の更新においては、
【数１１１】
0005842992.tif 000217
を用い、前記各粒子に係る力は、上記（式１２）の右辺第２項で示される圧力勾配項、右辺第３項で示される粘性応力項、右辺第４項で示される重力項、右辺第５項は粒子間距離が近づきすぎないようにするための力であり、右辺第６項は表面張力であり、ここで、
【数１１２】
0005842992.tif 000218
【数１１３】
0005842992.tif 000219
【数１１４】
0005842992.tif 000220
また、下付き添え字は粒子の番号を表し、
0005842992.tif 000221
：i番目の粒子の位置ベクトル、速度ベクトル、密度、圧力であり、
0005842992.tif 000222
：j番目の質量であり、
0005842992.tif 000223
：粘性項を計算するために導入した定数であり、
【数１１５】
0005842992.tif 000224
であり
0005842992.tif 000225
【数１１６】
0005842992.tif 000226
ｙ＝０の平面が固体の壁面の場合、
0005842992.tif 000227
であり、
ｇ：重力加速度
ｙ：ｙ方向の単位ベクトル
ｈ：カーネル関数の値が非ゼロとなる球形領域の半径である影響半径
Ｗ：ＳＰＨ法で用いられるカーネル関数の重ね合わせ
0005842992.tif 000228
は
【数１１７】
0005842992.tif 000229
の再標準化行列であり、
【数１１８】
0005842992.tif 000230
とおくと、
0005842992.tif 000231
0005842992.tif 000232
であり、
0005842992.tif 000233
は粒子ｉ，ｊ間での１次元リーマン問題を解いて求めた時空間の中間値であり、以下のよ
うに決定し、
まず、粒子ｉ，ｊに関して、以下の特性量を定め、
【数１１９】
0005842992.tif 000234
【数１２０】
0005842992.tif 000235
【数１２１】
0005842992.tif 000236
【数１２２】
0005842992.tif 000237
さらに、各勾配を以下のよう計算し、
【数１２３】
0005842992.tif 000238
【数１２４】
0005842992.tif 000239
【数１２５】
0005842992.tif 000240
【数１２６】
0005842992.tif 000241
【数１２７】
0005842992.tif 000242
【数１２８】
0005842992.tif 000243
ここで、
0005842992.tif 000244
0005842992.tif 000245
とし、上記の量を用いて、
0005842992.tif 000246
は以下のように決定する
【数１２９】
0005842992.tif 000247
【数１３０】
0005842992.tif 000248
【数１３１】
0005842992.tif 000249
【数１３２】
0005842992.tif 000250
【数１３３】
0005842992.tif 000251
ことを特徴とするシミュレーション装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、シミュレーションプログラム、シミュレーション方法及びシミュレーション装置に関する。
【背景技術】
【０００２】
近年の計算機の演算能力の向上につれて、計算機シミュレーションの手法も発展してきた。その結果として、様々な応用分野に計算機シミュレーションが用いられるようになってきた。【０００３】
流体や弾性体等の連続体の問題を解く数値計算手法としては、格子をベースにして微分方程式の近似解を求解する有限差分法、有限要素法、又は有限体積法等が多く用いられてきた。近年では数値計算をＣＡＥ（Computer Aided Engineering）などの応用分野で活用するため、これらの数値計算手法も発展し、流体と構造物が相互作用する問題が解かれるようになってきた。【０００４】
しかしながら、数値計算に格子を用いる手法では、自由表面などの界面の存在する問題や、流体・構造連成問題などの移動境界が発生する場合には、取り扱いが複雑なため、プログラム作成が困難である場合が多い。【０００５】
これに対して数値計算に格子を用いない粒子法（ＭＰＳ（Moving Particle Semi-implicit）法、ＳＰＨ（Smoothed Particle Hydrodynamics）法等）では、移動境界の取り扱いに特別な処置を必要としない。それゆえ、近年粒子法は広く用いられるようになっている。【０００６】
図１は、粒子法の問題点を説明するための図である。
自由表面などの、変形し移動する境界を容易に扱いたいがために開発された粒子法であるが、離散化方法から、連続体が、粒子１１から構成されるただの粒子群となってしまっているため、図１に示すように、連続体の境界面がどこかはっきりしなくなる。そのため、粒子法では表面張力などの境界を陽に扱う必要がある問題に対しては、統一的手法は出来ていない。【０００７】
差分法や有限要素法では、ＣＳＦ（Continuum Surface Force）モデルと呼ばれる表面張力の計算方法がある（例えば、非特許文献１参照。）。ＣＳＦモデルの粒子法への応用がX.Y. Hu, N.A. Adamsにより行われ、粒子法による表面張力の計算が可能となった（例えば、非特許文献２参照。）。また、他にも粒子間の相互作用として、表面張力を決定する手法（例えば、特許文献１参照。）や、近傍粒子の重心位置からのずれから表面張力を求める手法がある（例えば、非特許文献３、４参照。）。【０００８】
ここで、粒子法に適用されているＣＳＦモデルによる計算方法の概略を説明する。
図２は、影響半径、およびカーネル関数の重ねあわせによる連続関数を示す図である。
まず、表面張力の算出に用いるcolor関数＜Ｃi＞を以下のように定義する。
【数１】
0005842992.tif 000002
【０００９】
ここで、下付き添え字iはi番目の粒子１１の値を表わし、ｊはｊ番目の粒子１１の値を表わす。
0005842992.tif 000003
はそれぞれ粒子１１の位置ベクトル、粒子１１の質量および密度であり、各粒子１１に値が与えられている。
0005842992.tif 000004
はＳＰＨ法で用いられるカーネル関数２２の重ね合わせで、３次のスプライン関数２３などがよく用いられる。
0005842992.tif 000005
はカーネル関数２２の値が非ゼロとなる球形領域の半径を表し、図２に示したように影響半径２１と呼ばれる。【００１０】
color関数の勾配（▽Ｃi）はＳＰＨ法の標準的な手法であり、以下のように計算される。
【数２】
0005842992.tif 000006
【００１１】
上記（式２）を用いて粒子iにかかる表面張力を以下のように定義する。
【数３】
0005842992.tif 000007
【００１２】
ここで、
【数４】
0005842992.tif 000008
である。
0005842992.tif 000009
は表面張力係数、
0005842992.tif 000010
は次元、
0005842992.tif 000011
は単位行列である。【００１３】
ＣＳＦモデルによる計算方法は、上記（式３）により得られた表面張力を用いて時間発展させる計算方法である。【００１４】
図３乃至図７は、従来のＣＳＦモデルによる計算方法を用いた結果の経時的変化を示す図である。
図３乃至図７に示すように、ＣＳＦモデルによる計算方法は、流体の塊が表面張力の影響を受けながら重力により落下し、平面に衝突する計算方法である。図３乃至図７に示した挙動の計算は、非特許文献２に開示されている計算方法を用いた。図３乃至図７においては、ｙ＝０が壁面を表すように鉛直方向をｙ軸に設定しており、ｙ＜０となる位置には流体３１が侵入できないようになっている。【００１５】
図３に示すように、流体３１の初期形状は、０．０５[m]×０．０５[m]の正方形である。また、
【数５】
0005842992.tif 000012
[N/m2]、
【数６】
0005842992.tif 000013
[kg/m3]、
【数７】
0005842992.tif 000014
[m]、
【数８】
0005842992.tif 000015
[m/s]、
【数９】
0005842992.tif 000016
[s]、粘性率
0005842992.tif 000017
は０．００１[m2/s]、重力加速度
【数１０】
0005842992.tif 000018
[m/s2]、カーネル関数は３次のスプライン関数を用いた。初期状態の粒子位置は均等格子状におき（一辺当たり３３個で総数は１０８９個)、速度はすべての粒子で０を与える。【００１６】
計算結果を見ると、図４に示すように、流体３１は時刻０．０４[s]では滑らかな円形分布になり、図５に示すように、時刻０．１４[s]で壁と衝突し大きく変形する。
【先行技術文献】
【特許文献】
【００１７】
【特許文献１】特開２００８−１１１６７５号公報
【非特許文献】
【００１８】
【非特許文献１】B. Lafaurie, C. Nardone, R. Scardovelli, S. Zaleski, and G. Zanetti, “Modelling Merging and Fragmentation in Multiphase Flows with SuRFER”, Journal of Computaional Physics, 113, 134, 147, (1994)
【非特許文献２】X.Y. Hu, N.A. Adams, “A multi-phase SPH method for macroscopic and mesoscopic flows”, Journal of Computaional Physics, 213, pp. 844-861 (2006)
【非特許文献３】T. Hongo, M. Shigeta, S.  Izawa, and Y. Fukunishi, “３次元非圧縮ＳＰＨ法における気液界面に作用する表面張力モデル”, 第２３回数値流体力学シンポジウム, A8-5 (2009)
【非特許文献４】M. Agawa, M. Shigeta, S.  Izawa, and Y. Fukunishi, “斜面上を流下する液体の非圧縮ＳＰＨシミュレーション”, 第２３回数値流体力学シンポジウム, A9-4 (2009)
【発明の開示】
【００１９】
しかしながら、上述のような従来の表面張力計算手法においては、衝突した流体３１が大きく変形した場合に、流体３１において不自然な速度分布が発生するなどの挙動を示す場合がある、という問題点がある。
【００２０】
上述の例では、図６に示すように、時刻０．７[s]では形状が崩れ、図７に示すように、不自然な分布の速度分布３２が発生している。
【００２１】
１つの側面では、本発明は、不自然な分布の速度分布が発生するなどの挙動を抑制したシミュレーション結果を出力することが可能なシミュレーションプログラム、シミュレーション方法及びシミュレーション装置を提供することを目的とする。
【００２２】
１つの案では、流体の経時的な形状変化をシミュレーションする場合に、不自然な分布の速度分布が発生する挙動を抑制するシミュレーションプログラムにおいて、
コンピュータに、
格納部から、前記流体を形成する粒子の位置ベクトル及び速度ベクトルを取得し、
前記流体を粒子の集合で表現した流体モデルについて、取得した前記位置ベクトル及び速度ベクトルを用いて、前記流体を形成する各粒子に係る力および該各粒子に係る表面張力を算出し、
前記速度ベクトルに算出した各粒子に係る力及び各粒子に係る表面張力を加えて該各粒子の速度を更新し、
更新した該各粒子の速度に基づいて、前記流体の形状変化を算出し、
算出された前記形状変化を出力部に出力する
処理を実行させるシミュレーションプログラムであって、
前記該各粒子の速度の更新においては、
0005842992.tif 000019
を用い、前記各粒子に係る力は、上記（式１２）の右辺第２項で示される圧力勾配項、右辺第３項で示される粘性応力項、右辺第４項で示される重力項、右辺第５項は粒子間距離が近づきすぎないようにするための力であり、右辺第６項は表面張力であり、ここで、
0005842992.tif 000020
0005842992.tif 000021
0005842992.tif 000022
また、下付き添え字は粒子の番号を表し、
0005842992.tif 000023
：i番目の粒子の位置ベクトル、速度ベクトル、密度、圧力であり、
0005842992.tif 000024
：j番目の質量であり、
0005842992.tif 000025
：粘性項を計算するために導入した定数であり、
0005842992.tif 000026
であり
0005842992.tif 000027
0005842992.tif 000028
ｙ＝０の平面が固体の壁面の場合、
0005842992.tif 000029
であり、
ｇ：重力加速度
ｙ：ｙ方向の単位ベクトル
ｈ：カーネル関数の値が非ゼロとなる球形領域の半径である影響半径
Ｗ：ＳＰＨ法で用いられるカーネル関数の重ね合わせ
0005842992.tif 000030
は
0005842992.tif 000031
の再標準化行列であり、
0005842992.tif 000032
とおくと、
0005842992.tif 000033
0005842992.tif 000034
であり、
0005842992.tif 000035
は粒子ｉ，ｊ間での１次元リーマン問題を解いて求めた時空間の中間値であり、以下のように決定し、
まず、粒子ｉ，ｊに関して、以下の特性量を定め、
0005842992.tif 000036
0005842992.tif 000037
0005842992.tif 000038
0005842992.tif 000039
さらに、各勾配を以下のよう計算し、
0005842992.tif 000040
0005842992.tif 000041
0005842992.tif 000042
0005842992.tif 000043
0005842992.tif 000044
0005842992.tif 000045
ここで、
0005842992.tif 000046
0005842992.tif 000047
とし、上記の量を用いて、
0005842992.tif 000048
は以下のように決定する
0005842992.tif 000049
0005842992.tif 000050
0005842992.tif 000051
0005842992.tif 000052
0005842992.tif 000053
というシミュレーションプログラムが提供される。
【００２３】
不自然な分布の速度分布が発生するなどの挙動を抑制したシミュレーション結果を出力することができる。
【図面の簡単な説明】
【００２４】
【図１】粒子法の問題点を説明するための図である。
【図２】影響半径、およびカーネル関数の重ねあわせによる連続関数を示す図である。
【図３】従来のＣＳＦモデルによる計算方法を用いた結果の経時的変化（時刻＝０[s]）を示す図である。
【図４】従来のＣＳＦモデルによる計算方法を用いた結果の経時的変化（時刻＝０．０４[s]）を示す図である。
【図５】従来のＣＳＦモデルによる計算方法を用いた結果の経時的変化（時刻＝０．１４[s]）を示す図である。
【図６】従来のＣＳＦモデルによる計算方法を用いた結果の経時的変化（時刻＝０．７[s]）を示す図である。
【図７】従来のＣＳＦモデルによる計算方法を用いた結果の経時的変化（時刻＝０．７[s]：不自然な速度分布発生）を示す図である。
【図８】表面張力係数を説明するための図である。
【図９】本発明を適用したシミュレーション装置の構成例を示す図である。
【図１０】表面エネルギーΧ（カイ）を幅εで滑らかに近似した関数を説明するための図である。
【図１１】２次元流体の例を示す図である。
【図１２】粒子が壁面に衝突した状態を説明するための図である。
【図１３】本発明を適用したシミュレーション処理の流れを示すフローチャートである。
【図１４】流体の形状について、本発明を適用したシミュレーション処理を用いた結果の経時的変化（時刻＝０[s]）を示す図である。
【図１５】流体の形状について、本発明を適用したシミュレーション処理を用いた結果の経時的変化（時刻＝０．０４[s]）を示す図である。
【図１６】流体の形状について、本発明を適用したシミュレーション処理を用いた結果の経時的変化（時刻＝０．１４[s]）を示す図である。
【図１７】流体の形状について、本発明を適用したシミュレーション処理を用いた結果の経時的変化（時刻＝０．７[s]）を示す図である。
【図１８】流体の形状について、本発明を適用したシミュレーション処理を用いた結果の経時的変化（時刻＝０．７[s]：不自然な速度分布の発生無し）を示す図である。
【図１９】第２の実施の形態における近似関数の効果を示す図である。
【図２０】情報処理装置の構成図である。
【発明を実施するための形態】
【００２５】
以下、本発明の実施の形態について、図面を参照しながら詳細に説明する。
本発明を適用した実施の形態は、コンピュータに実行させるシミュレーションプログラム、シミュレーション方法及びシミュレーション装置であって、シミュレーションの対象である流体を粒子の集合として捉え、変分を基にした粒子間の相互作用力である表面張力を、所定の近似関数を用いて算出し、その算出した表面張力を用いることにより、流体の形状変化をシミュレーションする。
【００２６】
（第１の実施の形態）
図８は、表面張力係数を説明するための図である。
一般に流体の表面エネルギーは以下のように定義できる。
【数１１】
0005842992.tif 000054
【００２７】
そして、上記の表面エネルギーを粒子法により離散化して計算し、第一変分より流体粒子にかかる表面張力を求めることができる。
【００２８】
ここで、上記（式４）の右辺の積分領域は、流体の境界面全域であり、
0005842992.tif 000055
は流体が空気に触れている部分では１、それ以外の部分では０の値をとる。
0005842992.tif 000056
はそれぞれ、図８に示すように、空気に触れている部分の表面張力係数、固体に触れている部分の表面張力係数である。
【００２９】
図９は、本発明を適用したシミュレーション装置の構成例を示す図である。
図９において、シミュレーション装置９００は、処理部９０１、格納部９０２及び出力部９０３を備える。そして、シミュレーション装置９００は、後述する近似間数を用い、流体等の連続体を粒子の集まりとして表現する粒子法による数値計算を行って表面張力を算出する。その計算結果を用いて連続体の形状変化を算出し、シミュレーション結果９１２を出力する。
【００３０】
格納部９０２は、本発明を適用したシミュレーションプログラムを実行する場合に用いられる各計算式の情報を格納する。
【００３１】
処理部９０１は、第１および第２の実施の形態として後述する、本発明を適用したシミュレーション処理を実行する。
出力部９０３は、処理部９０１が実行したシミュレーション結果９１２を出力する。
【００３２】
まず、近似関数を用いて表面張力を算出する手法について説明する。
上述の（式４）の積分を、下記（式５）の近似関数のように定義する。
【数１２】
0005842992.tif 000057
【００３３】
ここで、
0005842992.tif 000058
は定数であり、表面張力波の振動周期と合うように決められる。
0005842992.tif 000059
は空間の次元、
0005842992.tif 000060
はそれぞれ、液体が空気と触れているときの表面張力と液体が固体と触れているときの表面張力である。
0005842992.tif 000061
は粒子の位置から決まる関数で、
【数１３】
0005842992.tif 000062
【数１４】
0005842992.tif 000063
【数１５】
0005842992.tif 000064
とする。
【００３４】
図１０は、表面エネルギーΧ（カイ）を幅εで滑らかに近似した関数を説明するための図である。
【００３５】
（式５）中の
0005842992.tif 000065
は（式４）の
0005842992.tif 000066
を幅
0005842992.tif 000067
で、図１０に示したように滑らかに近似した関数であり、
0005842992.tif 000068
は粒子iと固体境界との距離である。
0005842992.tif 000069
が面積の次元をもち、さらにＳＰＨ法では
0005842992.tif 000070
は流体の内部の粒子では１に近く、表面近くになると値が小さくなり、影響半径ｈ内に１つも粒子がなくなると、
0005842992.tif 000071
という値となる。
【００３６】
（式５）のように定義した近似関数は、流体表面の粒子の長さと同じ次元を持つ量を足し合わせたものであり、（式４）の表面積分を近似すると見なす。
【００３７】
（式５）で粒子法により離散化された表面エネルギーの勾配を計算し、各粒子にかかる正味の表面張力を導出すると以下のようになる。
【数１６】
0005842992.tif 000072
以上、近似関数を用いた表面張力の算出について説明した。
【００３８】
次に、上述のようにして算出した表面張力を用いて、流体の形状変化の算出について説明する。
【００３９】
図１１は、２次元流体の例を示す図であり、壁面となる平面から鉛直方向にｙ軸を設定するとともに、壁面に沿ってｘ軸を設定している。
【００４０】
図１１のような状況下、すなわちｙ＝０の固体１１１の平面の上に置かれた流体１１２での、非圧縮性粘性流体の運動方程式を考える。
【数１７】
0005842992.tif 000073
【数１８】
0005842992.tif 000074
【数１９】
0005842992.tif 000075
において、（式７）は質量保存則、（式８）は運動量保存則、（式９）は状態方程式である。
【数２０】
0005842992.tif 000076
はそれぞれ、流体１１２の密度場、速度場、圧力場、音速である。
0005842992.tif 000077
は粘性応力テンソルで、流体１１２の粘性係数を
0005842992.tif 000078
（定数）とすると、
【数２１】
0005842992.tif 000079
である。
【００４１】
（式７）は、流体１１２が集まってくるような速度場の場合に密度が上昇する効果を表わし、逆に流体１１２が離れていくような速度場の場合に密度が低下する効果を表わす。（式８）の右辺第１項は圧力勾配項で、流体１１２が圧力の大きい部分から圧力の小さい部分へ向かって力が発生する効果を表す。右辺第２項は粘性応力項で、流体の流れにブレーキがかかるような効果をあらわす。右辺第３項は重力項である。
0005842992.tif 000080
は重力加速度、
0005842992.tif 000081
はy方向の単位ベクトルである。
【００４２】
本第１の実施の形態では（式９）のような密度と圧力のみの関係を用いているが、一般の温度や内部エネルギー、エントロピー等を用いた状態方程式を用いても良い。
【００４３】
また、空気と接触部には
0005842992.tif 000082
の表面張力係数（定数）、固体１１１との接触部には
0005842992.tif 000083
の表面張力係数（定数）がかかっているものとする。
【００４４】
そして、（式７）乃至（式９）を、ＳＰＨ法を用いて以下のように離散化し、上述した表面張力項を導入する。
【数２２】
0005842992.tif 000084
【００４５】
ここで、（式１０）で計算した
0005842992.tif 000085
のｘ成分、ｙ成分をそれぞれ
0005842992.tif 000086
，
0005842992.tif 000087
とし、i番目の速度の成分を
0005842992.tif 000088
として、以下の計算を行う。
【数２３】
0005842992.tif 000089
【数２４】
0005842992.tif 000090
【数２５】
0005842992.tif 000091
【数２６】
0005842992.tif 000092
【数２７】
0005842992.tif 000093
【００４６】
ここで、下付き添え字は粒子の番号を表す。すなわちi番目の粒子の位置ベクトル、速度ベクトル、密度、圧力はそれぞれ
0005842992.tif 000094
である。
0005842992.tif 000095
はj番目の質量である。
0005842992.tif 000096
は粘性項を計算するために導入した定数である。
【００４７】
（式１２）は（式８）の運動方程式を粒子法により離散化したもので、右辺第２項は圧力勾配項、第３項は粘性応力項を表したものである。（式１２）の右辺第５項は粒子間距離が近づきすぎないようにするための力で、
0005842992.tif 000097
を定数とすれば、
【数２８】
0005842992.tif 000098
となる。
【００４８】
右辺第６項
0005842992.tif 000099
は、上述したような表面張力による力で、（式６）が適用されるが、上述のシミュレーションの場合ｙ＝０の平面が固体１１１の壁面であるため、
【数２９】
0005842992.tif 000100
である。
【００４９】
図１２は、粒子が壁面に衝突した状態を説明するための図である。
上述の（式１２）および（式１４）は、粒子の位置を移動させ、ｙ＝０の平面の位置にある壁面を通過しようとする場合に、図１２の上に示したように壁面をそのまま通過するのではなく、図１２の下に示したように壁面に衝突した瞬間に止まるという現実に即した修正を表す。
【００５０】
ここで、
0005842992.tif 000101
は
【数３０】
0005842992.tif 000102
の再標準化行列である。さらに、
【数３１】
0005842992.tif 000103
とおく。
【数３２】
0005842992.tif 000104
，
【数３３】
0005842992.tif 000105
である。
0005842992.tif 000106
は粒子ｉ，ｊ間での１次元リーマン問題を解いて求めた時空間の中間値である。具体的には以下のように決定する。
【００５１】
まず、粒子ｉ，ｊに関して、以下の特性量を定める。
【数３４】
0005842992.tif 000107
【数３５】
0005842992.tif 000108
【数３６】
0005842992.tif 000109
【数３７】
0005842992.tif 000110
【００５２】
さらに、各勾配を以下のよう計算する。
【数３８】
0005842992.tif 000111
【数３９】
0005842992.tif 000112
【数４０】
0005842992.tif 000113
【数４１】
0005842992.tif 000114
【数４２】
0005842992.tif 000115
【数４３】
0005842992.tif 000116
【００５３】
ここで、
【数４４】
0005842992.tif 000117
【数４５】
0005842992.tif 000118
とし、上記の量を用いて、
0005842992.tif 000119
は以下のように決定する。
【数４６】
0005842992.tif 000120
【数４７】
0005842992.tif 000121
【数４８】
0005842992.tif 000122
【数４９】
0005842992.tif 000123
【数５０】
0005842992.tif 000124
【００５４】
図１３は、本発明を適用したシミュレーション処理の流れを示すフローチャートである。
【００５５】
まず、ステップＳ１３０１において、シミュレーションの対象とする流体の物理量のデータとして、各粒子の位置ベクトル、速度ベクトルを取得し、ステップＳ１３０２において、上述の（式９）を用いて、現在の速度を用いて粒子の位置を時間ｄｔ／２更新する。
【００５６】
そして、（式１１）におけるｙの値による場合分けにおいて、ｙ＜０となる結果、粒子が壁面を通り抜けると判断された場合（ステップＳ１３０３：Ｙｅｓ）は、ステップＳ１３０４において、上述の（式１１）を用いて、粒子が壁面を通り抜けないよう、すなわち、ｙ≧０となるように粒子の位置と速度を補正する。
【００５７】
一方、（式１１）におけるｙの値による場合分けにおいて、ｙ≧０となる結果、粒子が壁面を通り抜けていないと判断された場合（ステップＳ１３０３：Ｎｏ）は、粒子の位置と速度は補正しない。
【００５８】
次に、ステップＳ１３０５において、上述の（式１２）右辺の第２項乃至第５項を算出することにより、物理モデルにより流体粒子にかかる力を算出する。
【００５９】
そして、ステップＳ１３０６において、上述の（式１２）右辺の第６項を算出することにより、各流体粒子に係る表面張力を算出する。
【００６０】
各境界粒子に対する表面張力の計算が終了したら、ステップＳ１３０７において、流体粒子の速度を更新する。
【００６１】
次に、ステップＳ１３０８において、上述の（式１３）を用いて、ステップＳ１３０７で更新した速度を用いて流体粒子の位置を時間ｄｔ／２経過後の位置に更新する。
【００６２】
そして、（式１４）におけるｙの値による場合分けにおいて、ｙ＜０となる結果、粒子が壁面を通り抜けると判断された場合（ステップＳ１３０９：Ｙｅｓ）は、ステップＳ１３１０において、上述の（式１４）を用いて、粒子が壁面を通り抜けないよう、すなわち、ｙ≧０となるように粒子の位置と速度を補正する。
【００６３】
一方、（式１４）におけるｙの値による場合分けにおいて、ｙ≧０となる結果、粒子が壁面を通り抜けていないと判断された場合（ステップＳ１３０３：Ｎｏ）は、粒子の位置と速度は補正しない。
【００６４】
その後、ステップＳ１３１１において、流体粒子の速度を更新し、上記ステップＳ１３０２乃至ステップＳ１３１１を繰り返し実行する。
【００６５】
以上、説明したシミュレーション処理を実行すると、流体の形状は図１４乃至図１８に示したようになる。
【００６６】
図１４乃至図１８は、流体の形状について、本発明を適用したシミュレーション処理を用いた結果の経時的変化を示す図である。なお、図３乃至図７の従来技術での計算と同様の計算条件である。
【００６７】
図１４に示すように、流体１４１の初期形状は、０．０５[m]×０．０５[m]の正方形である。また、それぞれの初期条件を、
【数５１】
0005842992.tif 000125
[N/m2]、
【数５２】
0005842992.tif 000126
[kg/m3]、
【数５３】
0005842992.tif 000127
[m]、
【数５４】
0005842992.tif 000128
[m/s]、
【数５５】
0005842992.tif 000129
[s]、
【数５６】
0005842992.tif 000130
、粘性率
0005842992.tif 000131
は０．００１[m2/s]、重力加速度
【数５７】
0005842992.tif 000132
[m/s2]、
【数５８】
0005842992.tif 000133
、
【数５９】
0005842992.tif 000134
[m]とした。
【００６８】
初期状態の粒子は均等格子状に配置し（一辺当たり３３個で総数は１０８９個)、速度はすべての粒子で０と密度は１０００[kg/m3]を与える。
【００６９】
計算結果を見ると、図１５に示すように、流体１４１は時刻０．０４[s]では図４と同様に均一な円形分布になる。そして、図１６に示すように、時刻０．１４[s]で壁と衝突し大きく変形する。
【００７０】
図１７に示したように、時刻０．７[s]では、図６に示した従来手法の計算結果とは異なり、流体１４１の形状は崩れず、図１８に示すように、速度分布も落ち着き定常状態が維持されている。
【００７１】
（第２の実施の形態）
本発明を適用した第２の実施の形態では、上述の第１の実施の形態で、（式４）の積分を（式５）の近似関数のように近似した代わりに、（式４）の積分を下記（式３１）の近似関数のように近似する。第１の実施の形態中の圧力「ｐ１」を「ｐ２」に置き換えることにより、第１の実施の形態中の説明を第２の実施の形態に適用できる。
【数６０】
0005842992.tif 000135
【００７２】
ここで
【数６１】
0005842992.tif 000136
【数６２】
0005842992.tif 000137
【数６３】
0005842992.tif 000138
である。
0005842992.tif 000139
も面積と同様の次元を持つため、（式６）は表面エネルギーを近似したものと見なす。
0005842992.tif 000140
は再標準化行列であり、粒子の非一様な分布における勾配計算を補正する。
【００７３】
図１９は、第２の実施の形態における近似関数の効果を示す図である。
上記（式３１）中の
0005842992.tif 000141
１９２は、
0005842992.tif 000142
１９１を原点付近で滑らかにした近似関数で、
0005842992.tif 000143
極小値の大きさが、
0005842992.tif 000144
とほぼ同じであり、「０．１」となっている。
【００７４】
（式６）のエネルギーを基にして粒子ｋにかかる正味の力を導出すると以下のようになる。
【数６４】
0005842992.tif 000145
0005842992.tif 000146
はそれぞれ、ｘ，ｙ，ｚ方向の単位ベクトルである。
【００７５】
上述したように、本発明は粒子法の計算で、表面張力が効果的な流体運動の計算に対して適用可能である。特に溶けた金属の流し込みのシミュレーションや、樹脂の流し込みのシミュレーション等には効果的である。また、壁面等の固体に衝突するような流体の形状変化をシミュレーションする際に、従来のように壁面等に粒子を置く設定を行う必要が無くなる。
【００７６】
図９のシミュレーション装置は、例えば、図２０に示すような情報処理装置（コンピュータ）を用いて実現することが可能である。図２０の情報処理装置は、ＣＰＵ（Central Processing Unit）２００１、メモリ２００２、入力装置２００３、出力装置２００４、外部記憶装置２００５、媒体駆動装置２００６及びネットワーク接続装置２００７を備える。これらはバス２００８により互いに接続されている。
【００７７】
メモリ２００２は、例えば、ＲＯＭ（Read Only Memory）、ＲＡＭ（Random Access Memory）、フラッシュメモリ等の記憶装置であり、シミュレーション処理に用いられるプログラム及びデータを格納する。例えば、ＣＰＵ２００１は、メモリ２００２を利用してプログラムを実行することにより、上述のシミュレーション処理を行う。メモリ２００２は、図９の格納部９０２としても使用できる。
【００７８】
入力装置２００３は、例えば、キーボード、ポインティングデバイス等であり、オペレータからの指示や情報の入力に用いられる。出力装置２００４は、例えば、表示装置、プリンタ、スピーカ等であり、オペレータへの問い合わせや処理結果の出力に用いられる。出力装置２００４は、図９の出力部９０３としても使用できる。
【００７９】
外部記憶装置２００５は、例えば、磁気ディスク装置、光ディスク装置、光磁気ディスク装置、テープ装置等である。この外部記憶装置２００５には、ハードディスクドライブも含まれる。情報処理装置は、この外部記憶装置２００５にプログラム及びデータを格納しておき、それらをメモリ２００２にロードして使用することができる。
【００８０】
媒体駆動装置２００６は、可搬型記録媒体２００９を駆動し、その記録内容にアクセスする。可搬型記録媒体２００９は、メモリデバイス、フレキシブルディスク、光ディスク、光磁気ディスク等である。この可搬型記録媒体２００９には、Compact Disk Read Only Memory（ＣＤ−ＲＯＭ）、Digital Versatile Disk（ＤＶＤ）、Universal Serial Bus（ＵＳＢ）メモリ等も含まれる。オペレータは、この可搬型記録媒体２００９にプログラム及びデータを格納しておき、それらをメモリ２００２にロードして使用することができる。
【００８１】
このように、シミュレーション処理に用いられるプログラム及びデータを格納するコンピュータ読み取り可能な記録媒体には、メモリ２００２、外部記憶装置２００５、及び可搬型記録媒体２００９のような、物理的な（非一時的な）記録媒体が含まれる。
【００８２】
ネットワーク接続装置２００７は、通信ネットワーク２０１０に接続され、通信に伴うデータ変換を行う通信インタフェースである。情報処理装置は、プログラム及びデータを外部の装置からネットワーク接続装置２００７を介して受け取り、それらをメモリ２００２にロードして使用することができる。ネットワーク接続装置２００７は、図９の出力部９０３としても使用できる。
【００８３】
開示の実施形態とその利点について詳しく説明したが、当業者は、特許請求の範囲に明確に記載した本発明の範囲から逸脱することなく、様々な変更、追加、省略をすることができるであろう。
【００８４】
例えば、（式５）または（式３１）のように定義した近似関数の代わりに、適切な近似関数を用いることもできる。
【００８５】
すなわち、本発明は、以上に述べた実施の形態に限定されるものではなく、本発明の要旨を逸脱しない範囲内で種々の構成または形状を取ることができる。
【００８６】
そして、開示されたシミュレーションプログラム、シミュレーション方法及びシミュレーション装置によれば、シミュレーション対象である流体について、不自然な速度分布が発生するなどの挙動を示すことなく、適切なシミュレーション結果を出力する、という効果を奏する。
【図１】
0005842992.tif 000252
【図２】
0005842992.tif 000253
【図３】
0005842992.tif 000254
【図４】
0005842992.tif 000255
【図５】
0005842992.tif 000256
【図６】
0005842992.tif 000257
【図７】
0005842992.tif 000258
【図８】
0005842992.tif 000259
【図９】
0005842992.tif 000260
【図１０】
0005842992.tif 000261
【図１１】
0005842992.tif 000262
【図１２】
0005842992.tif 000263
【図１３】
0005842992.tif 000264
【図１４】
0005842992.tif 000265
【図１５】
0005842992.tif 000266
【図１６】
0005842992.tif 000267
【図１７】
0005842992.tif 000268
【図１８】
0005842992.tif 000269
【図１９】
0005842992.tif 000270
【図２０】
0005842992.tif 000271