(57)【要約】  粒子法を用いた運動シミュレーションにおいて、粒子群による大きなスケールの運動から各粒子の乱雑な運動が生成されることを防ぐ。格納部４０２は、現在時刻における複数の粒子の位置、速度、加速度、圧力、又は密度に付加された擾乱項を含み、前の時刻におけるそれらの粒子の位置、速度、及び密度（位置等）から現在時刻における位置等を求める計算式を格納する。処理部４０１は、上記擾乱項の値を求め、初期時刻における位置等を前の時刻における位置等として用い、擾乱項の値を上記計算式に適用することで、現在時刻における複数の粒子の位置等を求める。そして、得られた位置等を前の時刻における位置等として用いて複数の粒子の擾乱項の値と現在時刻における位置等を求める処理を繰り返す。出力部４０３は、初期時刻以降の複数の時刻における複数の粒子の位置をシミュレーション結果として出力する。
2012111082.tif 000002
【特許請求の範囲】
【請求項１】
格納部に格納された、現在時刻における複数の粒子の位置、速度、加速度、圧力、又は密度に付加された擾乱項の情報を用いて、該擾乱項の値を求め、
前記複数の粒子の前記擾乱項の値を前記格納部に格納し、
前記格納部に格納された、前記複数の粒子の前記擾乱項を含む計算式であって、前の時刻における該複数の粒子の位置、速度、及び密度から前記現在時刻における該複数の粒子の位置、速度、及び密度を求める該計算式の情報を用い、初期時刻における該複数の粒子の位置、速度、及び密度を該前の時刻における該複数の粒子の位置、速度、及び密度として用い、前記複数の粒子の前記擾乱項の値を該計算式に適用することで、該現在時刻における該複数の粒子の位置、速度、及び密度を求め、
得られた前記複数の粒子の位置、速度、及び密度を前記前の時刻における前記複数の粒子の位置、速度、及び密度として用いて該複数の粒子の前記擾乱項の値と前記現在時刻における該複数の粒子の位置、速度、及び密度を求める処理を繰り返すことで、前記初期時刻より後の複数の時刻における該複数の粒子の位置、速度、及び密度を順次求め、
前記複数の時刻における前記複数の粒子の位置をシミュレーション結果として出力する処理をコンピュータに実行させるプログラム。
【請求項２】
前記コンピュータは、前記現在時刻における各粒子の位置又は速度に付加された擾乱の向きを示す乱数を用いて、前記複数の粒子の位置又は速度に付加された前記擾乱項の値を求めることを特徴とする請求項１記載のプログラム。
【請求項３】
前記コンピュータは、前記現在時刻における各粒子の加速度、圧力、又は密度に付加された擾乱の大きさを示す乱数を用いて、前記複数の粒子の加速度、圧力、又は密度に付加された前記擾乱項の値を求めることを特徴とする請求項１記載のプログラム。
【請求項４】
前記コンピュータは、前記前の時刻における各粒子の位置、速度、及び密度を用いて該前の時刻における各粒子のエネルギーを求め、得られたエネルギーが前記初期時刻における該粒子のエネルギーより大きい場合、該粒子の擾乱項の値を零に設定することを特徴とする請求項１、２、又は３記載のプログラム。
【請求項５】
現在時刻における複数の粒子の位置、速度、加速度、圧力、又は密度に付加された擾乱項を含み、前の時刻における該複数の粒子の位置、速度、及び密度から該現在時刻における該複数の粒子の位置、速度、及び密度を求める計算式の情報を格納する格納部と、
前記複数の粒子の前記擾乱項の値を求め、該複数の粒子の該擾乱項の値を前記格納部に格納し、初期時刻における該複数の粒子の位置、速度、及び密度を前記前の時刻における前記複数の粒子の位置、速度、及び密度として用い、該複数の粒子の該擾乱項の値を前記計算式に適用することで、前記現在時刻における該複数の粒子の位置、速度、及び密度を求め、得られた前記複数の粒子の位置、速度、及び密度を前記前の時刻における前記複数の粒子の位置、速度、及び密度として用いて該複数の粒子の該擾乱項の値と該現在時刻における該複数の粒子の位置、速度、及び密度を求める処理を繰り返すことで、該初期時刻より後の複数の時刻における該複数の粒子の位置、速度、及び密度を順次求める処理部と、
前記複数の時刻における前記複数の粒子の位置をシミュレーション結果として出力する出力部と
を備えることを特徴とするシミュレーション装置。
【請求項６】
前記処理部は、前記現在時刻における各粒子の位置又は速度に付加された擾乱の向きを示す乱数を用いて、前記複数の粒子の位置又は速度に付加された前記擾乱項の値を求めることを特徴とする請求項５記載のシミュレーション装置。
【請求項７】
前記処理部は、前記現在時刻における各粒子の加速度、圧力、又は密度に付加された擾乱の大きさを示す乱数を用いて、前記複数の粒子の加速度、圧力、又は密度に付加された前記擾乱項の値を求めることを特徴とする請求項５記載のシミュレーション装置。
【請求項８】
前記処理部は、前記前の時刻における各粒子の位置、速度、及び密度を用いて該前の時刻における各粒子のエネルギーを求め、得られたエネルギーが前記初期時刻における該粒子のエネルギーより大きい場合、該粒子の擾乱項の値を零に設定することを特徴とする請求項５、６、又は７記載のシミュレーション装置。
【請求項９】
コンピュータにより実行されるシミュレーション方法であって、
格納部に格納された、現在時刻における複数の粒子の位置、速度、加速度、圧力、又は密度に付加された擾乱項の情報を用いて、該擾乱項の値を求め、
前記複数の粒子の前記擾乱項の値を前記格納部に格納し、
前記格納部に格納された、前記複数の粒子の前記擾乱項を含む計算式であって、前の時刻における該複数の粒子の位置、速度、及び密度から前記現在時刻における該複数の粒子の位置、速度、及び密度を求める該計算式の情報を用い、初期時刻における該複数の粒子の位置、速度、及び密度を該前の時刻における該複数の粒子の位置、速度、及び密度として用い、前記複数の粒子の前記擾乱項の値を該計算式に適用することで、該現在時刻における該複数の粒子の位置、速度、及び密度を求め、
得られた前記複数の粒子の位置、速度、及び密度を前記前の時刻における前記複数の粒子の位置、速度、及び密度として用いて該複数の粒子の前記擾乱項の値と前記現在時刻における該複数の粒子の位置、速度、及び密度を求める処理を繰り返すことで、前記初期時刻より後の複数の時刻における該複数の粒子の位置、速度、及び密度を順次求め、
前記複数の時刻における前記複数の粒子の位置をシミュレーション結果として出力することを特徴とするシミュレーション方法。
【請求項１０】
前記コンピュータは、前記現在時刻における各粒子の位置又は速度に付加された擾乱の向きを示す乱数を用いて、前記複数の粒子の位置又は速度に付加された前記擾乱項の値を求めることを特徴とする請求項９記載のシミュレーション方法。
【請求項１１】
前記コンピュータは、前記現在時刻における各粒子の加速度、圧力、又は密度に付加された擾乱の大きさを示す乱数を用いて、前記複数の粒子の加速度、圧力、又は密度に付加された前記擾乱項の値を求めることを特徴とする請求項９記載のシミュレーション方法。
【請求項１２】
前記コンピュータは、前記前の時刻における各粒子の位置、速度、及び密度を用いて該前の時刻における各粒子のエネルギーを求め、得られたエネルギーが前記初期時刻における該粒子のエネルギーより大きい場合、該粒子の擾乱項の値を零に設定することを特徴とする請求項９、１０、又は１１記載のシミュレーション方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、複数の時刻における複数の粒子の位置を求めるシミュレーション装置、シミュレーション方法、及びプログラムに関する。
【背景技術】
【０００２】
流体、弾性体等の連続体や粉粒体の運動を粒子群で表現するシミュレーション方法として粒子法が知られている。粒子法には、Moving Particle Semi-implicit （ＭＰＳ）法、Smoothed Particle Hydrodynamics （ＳＰＨ）法等が含まれる。この粒子法は、格子を用いる有限差分法や有限要素法等のシミュレーション方法に比べて、自由表面流等の移動境界を含む問題を扱いやすいため、近年、様々な分野で用いられるようになっている。【０００３】
粒子法は、自由表面流等の移動境界を含む問題を解くために開発された方法であるが、粒子の運動を安定に計算するために、人工粘性項と呼ばれる、粒子間の運動を抑える効果を持つ項が計算式に付加される。【０００４】
例えば、ＳＰＨ法では、運動量保存則と質量保存則が以下のように離散化される。【０００５】
【数１】
2012111082.tif 000003
【０００６】
式（１）及び（２）において、下付添え字ａはａ番目の粒子（粒子ａ）を表し、下付添え字ｂはｂ番目の粒子（粒子ｂ）を表す。
2012111082.tif 000004
ρa ，Ｐa ，ｍa は、それぞれ粒子ａの位置ベクトル、速度ベクトル、密度、圧力、及び質量を表す。粒子ｂについても同様の表記が用いられる。式（１）及び（２）の右辺の
2012111082.tif 000005
は、粒子ａ以外のすべての粒子についての総和を表す。【０００７】
Ｗはカーネル関数であり、粒子の分布から連続場を構成するのに用いられる。カーネル関数Ｗとしては、以下の３次のスプライン関数等がよく用いられる。【０００８】
【数２】
2012111082.tif 000006
【０００９】
ｈは粒子間の影響半径を表し、初期状態の平均粒子間隔の２倍から３倍程度の値がよく用いられる。φはカーネル関数Ｗの全空間積分量が１になるように調整された規格化因子であり、２次元の場合は０．７πｈ2 、３次元の場合はπｈ3 と決められる。【００１０】
式（１）の右辺第２項のΠabは粒子ａと粒子ｂの間に働く人工粘性項であり、例えば、以下のような形式が採用される。【００１１】
【数３】
2012111082.tif 000007
α、β、及びηは定数であり、ｃは音速を表す。【００１２】
式（１）は流体の運動量保存則を表す方程式であり、式（１）の右辺第１項は流体の圧力による力を表し、粒子同士が互いに押し合うような効果を持つ。人工粘性項Πabは、粒子同士が互いに近づこうとする際にブレーキとして運動量を減衰させ、粒子同士がすり抜けることを防ぐ効果を持つ。式（２）は流体の質量保存則を表す方程式であり、粒子同士が互いに近づく場合は密度が高くなり、遠ざかる場合は密度が小さくなる効果を持つ。【００１３】
しかしながら、人工粘性項を用いた式（１）−（３）の形式では、人工粘性項による運動の減衰の効果が目立ってしまう。そこで、人工粘性項による運動の減衰を極力抑えることを目的とした、リーマン問題の解を用いたGodunov SPH 法等も考案されている。【００１４】
また、流体と固体表面とのヌレ及びハジキ現象を再現する粒子型流体シミュレーション方法や、格子状に分割された空間に存在する粒子のデータに対する演算処理を効率良く実行するデータ処理方法も考案されている。
【先行技術文献】
【特許文献】
【００１５】
【特許文献１】特開平１０−１８５７５５号公報
【特許文献２】国際公開第２００８／０２０６３４号パンフレット
【非特許文献】
【００１６】
【非特許文献１】S. Koshizuka, H. Tamako and Y. Oka, “A Particle Method for Incompressible Viscous Flow with Fluid Fragmentation”, Computational Fluid Dynamics Journal, vol.4, no.1, p.29-46, April 1995
【非特許文献２】J.J. Monaghan, “Smoothed Particle Hydrodynamics”, Annual Review of Astronomy and Astrophysics, vol.30,p.543-574, 1992
【非特許文献３】S. Inutsuka, “Reformulation of Smoothed Particle Hydrodynamics with Riemann Solver”, Journal of Computational Physics 179, p.238-267, 2002
【発明の概要】
【発明が解決しようとする課題】
【００１７】
上述したGodunov SPH 法等を用いても、粒子群のような粒子直径の１０倍程度以上の規模の大きなスケールの運動から、各粒子の乱雑な運動が生成されてしまい、大きなスケールの整った流れが散逸してしまう結果、シミュレーションの精度が悪化してしまうという問題は、依然として残されている。【００１８】
例えば、水面波の伝播等の波動解析において、粒子法を用いて数値計算を進めて行くと、整った流れから各粒子の乱雑な運動が生成されてしまい、波高が現実よりも過剰に減衰してしまうことがある。図１は、このような整った流れを示している。【００１９】
図１の破線１０１〜１０３は、流体中の整った流れの流線を表し、ベクトル１２１〜１２５は、それぞれ粒子１１１〜１１５の速度ベクトルを表す。このとき、粒子１１１及び１１３は、粒子１１１に働くベクトル１３１及び粒子１１３に働くベクトル１３２が示すように、流線１０１〜１０３から外れる方向の力を粒子１１４から受けてしまう。このように、粒子直径の１０倍程度以上の規模の大きなスケールの整った流れからは、剛体球の玉突き衝突のような現象が起こり、各粒子の乱雑な運動が生成されるとともに、流線１０１〜１０３が示す整った運動が壊されてしまうものと考えられる。この問題は、系の運動エネルギーが保存しても起こる、粒子法に顕著な問題である。【００２０】
粒子法による水面波の散逸は、粒子直径の１０倍程度以上の規模の大きなスケールの整った流れから各粒子の乱雑な運動が生成されるプロセスと、生成された乱雑な運動が人工粘性項により散逸されるプロセスの２つのプロセスによるものと考えられる。Godunov SPH 法は、人工粘性項の効果を小さくするように開発された方法であるが、それでも小さなスケールでの乱雑な運動には散逸効果が現れ、運動を散逸させている。【００２１】
本発明の課題は、粒子法を用いた運動シミュレーションにおいて、粒子群による大きなスケールの運動から各粒子の乱雑な運動が生成されることを防ぐことである。
【課題を解決するための手段】
【００２２】
シミュレーション装置は、格納部、処理部、及び出力部を備える。
格納部は、現在時刻における複数の粒子の位置、速度、加速度、圧力、又は密度に付加された擾乱（じょうらん）項を含み、前の時刻におけるそれらの粒子の位置、速度、及び密度から現在時刻におけるそれらの粒子の位置、速度、及び密度を求める計算式の情報を格納する。【００２３】
処理部は、上記複数の粒子の擾乱項の値を求め、それらの粒子の擾乱項の値を格納部に格納する。次に、初期時刻におけるそれらの粒子の位置、速度、及び密度を前の時刻におけるそれらの粒子の位置、速度、及び密度として用い、それらの粒子の擾乱項の値を上記計算式に適用することで、現在時刻におけるそれらの粒子の位置、速度、及び密度を求める。【００２４】
そして、得られた上記複数の粒子の位置、速度、及び密度を前の時刻におけるそれらの粒子の位置、速度、及び密度として用いてそれらの粒子の擾乱項の値と現在時刻におけるそれらの粒子の位置、速度、及び密度を求める処理を繰り返す。これにより、初期時刻より後の複数の時刻におけるそれらの粒子の位置、速度、及び密度を順次求める。【００２５】
出力部は、上記複数の時刻における上記複数の粒子の位置をシミュレーション結果として出力する。
【発明の効果】
【００２６】
開示のシミュレーション装置によれば、粒子法を用いた運動シミュレーションにおいて、各粒子の位置、速度、加速度、圧力、又は密度に擾乱項を付加することで、粒子群による大きなスケールの運動から各粒子の乱雑な運動が生成されることが防止される。例えば、水面波等の問題に対するシミュレーションにおいて、整った流れから生成される乱雑な運動による波高の過剰な減衰を防止することができ、その結果、シミュレーション精度を向上させることが可能となる。
【図面の簡単な説明】
【００２７】
【図１】水面波シミュレーションにおける整った流れを示す図である。
【図２】粒子の速度に対する擾乱を示す図である。
【図３】粒子の位置に対する擾乱を示す図である。
【図４】シミュレーション装置の構成図である。
【図５】第１のシミュレーション処理のフローチャートである。
【図６】初期時刻における流体の２次元形状を示す図である。
【図７】擾乱項を付加しない場合のシミュレーション結果を示す図（その１）である。
【図８】擾乱項を付加しない場合のシミュレーション結果を示す図（その２）である。
【図９】擾乱項を付加した場合のシミュレーション結果を示す図（その１）である。
【図１０】擾乱項を付加した場合のシミュレーション結果を示す図（その２）である。
【図１１】第２のシミュレーション処理のフローチャートである。
【図１２】情報処理装置の構成図である。
【発明を実施するための形態】
【００２８】
以下、図面を参照しながら、実施形態を詳細に説明する。
粒子法の数値計算アルゴリズムには、各粒子の位置及び速度を更新する手順が含まれる。この手順としては、以下のような形式が採用されることが多い。
【００２９】
【数４】
2012111082.tif 000008
【００３０】
式（１１）及び（１２）において、下付添え字ａ及びｂはそれぞれ粒子ａ及びｂを表し、上付添え字ｎ及びｎ＋１はシミュレーションにおける時刻（時間ステップ）を表す。ｄｔは時間について離散化した際の時間刻み幅を表し、
2012111082.tif 000009
は、粒子ａが粒子ｂから受ける力を表す。
【００３１】
ここで、各粒子の速度に擾乱を与える場合、δを小さな実数又は乱数とし、Γを実数として、
2012111082.tif 000010
を単位ベクトルとして、式（１１）を以下のような形式に拡張する。
【００３２】
【数５】
2012111082.tif 000011
【００３３】
2012111082.tif 000012
の向きは、粒子群の速度分布や位置関係等、計算中の粒子の情報から確定的に決めても良く、乱数を用いて決めても良い。この場合、図２に示すように、各粒子の速度ベクトル２０１に擾乱ベクトル２０３が外力として加えられ、速度ベクトル２０２が生成される。
【００３４】
また、各粒子の位置に擾乱を与える場合、式（１２）を以下のような形式に拡張する。
【００３５】
【数６】
2012111082.tif 000013
【００３６】
この場合、図３に示すように、各粒子の位置ベクトルに擾乱ベクトル３０１が変位として加えられる。擾乱ベクトル３０１の向きは、計算中の粒子の情報から確定的に決めても良く、乱数を用いて決めても良い。
【００３７】
また、各粒子の加速度に擾乱を与える場合、式（１１）を以下のような形式に拡張する。
【００３８】
【数７】
2012111082.tif 000014
【００３９】
この場合、粒子ａが粒子ｂから受ける力に擾乱が加えられる。
さらに、流体計算を行なう際には各粒子の速度、位置、又は加速度に擾乱を与える代わりに、各粒子の圧力もしくは密度に擾乱を与えても良い。この場合、各粒子の擾乱を加えた圧力もしくは密度を以下のように決め、粒子aの圧力もしくは密度として用いる。
【００４０】
【数８】
2012111082.tif 000015
【００４１】
【数９】
2012111082.tif 000016
【００４２】
式（１６）について圧力の擾乱として付加されている実数δｐや式（１７）について密度の擾乱として付加されている実数δρは、粒子群の速度分布や位置関係等、計算中の粒子の情報から確定的に決めても良く、乱数を用いて決めても良い。
【００４３】
このように、各粒子の位置、速度、加速度、圧力、又は密度に擾乱を与えることで、人工粘性項により散逸してしまう各粒子の乱雑な運動を付加することができ、散逸された運動エネルギーを注入することができる。これにより、大きなスケールの運動から小さなスケールのランダムな運動へと運動エネルギーが流入してしまう現象を防ぐことが可能になる結果、シミュレーションの精度を向上させることができる。
【００４４】
粒子のランダムな運動では、人工粘性項による効果により運動エネルギーが散逸してしまうが、粒子のランダムな運動以外の大きなスケールの運動では、人工粘性項による効果は小さいと考えられる。さらに、微小なスケールのランダムな運動がある程度運動エネルギーを保有している場合は、粒子直径の１０倍程度以上の規模の大きなスケールの運動から運動エネルギーが流入しなくなると考えられる。したがって、各粒子の物理量に擾乱を与えて粒子のランダムな運動エネルギーを注入すれば、粒子直径の１０倍程度以上の規模の大きなスケールの運動を散逸させることなく、各粒子の位置を計算することができる。
【００４５】
図４は、実施形態のシミュレーション装置の構成例を示している。図４のシミュレーション装置は、処理部４０１、格納部４０２、及び出力部４０３を備え、初期条件４１１に基づいて粒子法による数値計算を行って、シミュレーション結果４１２を出力する。
【００４６】
格納部４０２は、現在時刻における複数の粒子についての物理量に付加された擾乱項を含み、前の時刻におけるそれらの粒子の位置、速度、及び密度から現在時刻におけるそれらの粒子の位置、速度、及び密度を求める計算式の情報を格納する。
なお、本実施形態では、各粒子の速度の物理量に対して擾乱を加えるシミュレーションについて説明をおこなうが、速度以外の物理量、例えば、位置、加速度、密度や圧力等に対して擾乱を加えてシミュレーションを行っても良い。
【００４７】
処理部４０１は、初期条件４１１に含まれる、初期時刻におけるそれぞれの粒子の位置、速度、及び密度を用いて、それぞれの粒子の速度についての擾乱項の値を求め、得られた擾乱項の値を格納部４０２に格納する。次に、初期時刻におけるそれぞれの粒子の位置、速度、及び密度を前の時刻におけるそれぞれの粒子の位置、速度、及び密度として用い、それぞれの粒子の速度についての擾乱項の値を計算式に適用することで、現在時刻におけるそれぞれの粒子の位置、速度、及び密度を求める。
【００４８】
そして、処理部４０１は、得られたそれぞれの粒子の位置、速度、及び密度を前の時刻におけるそれぞれの粒子の位置、速度、及び密度として用いて、それぞれの粒子の速度についての擾乱項の値と現在時刻におけるそれぞれの粒子の位置、速度、及び密度を求める処理を繰り返す。これにより、初期時刻より後の複数の時刻におけるそれぞれの粒子の位置、速度、及び密度が、順次求められる。
【００４９】
出力部４０３は、得られた複数の時刻におけるそれぞれの粒子の位置をシミュレーション結果４１２として出力する。
まず、以下の３次元流体の運動方程式を用いて、粒子法によりシミュレーションを行う場合を考える。
【００５０】
【数１０】
2012111082.tif 000017
【００５１】
式（２１）〜（２３）において、ρ（ｘ，ｙ，ｚ，ｔ）、ｐ（ｘ，ｙ，ｚ，ｔ）、及び
2012111082.tif 000018
は、それぞれ、流体の密度場、圧力場、及び速度場を現し、ｃは音速を表す。
2012111082.tif 000019
はベクトル微分演算子を表し、ｘ、ｙ、及びｚは空間座標を表し、ｔは時刻を表す。
2012111082.tif 000020
は、z方向（鉛直方向）の単位ベクトルを表し、ｇは重力加速度を表す。ρ0 は、圧力が０のときの密度を表す。
【００５２】
式（２１）は質量保存則を表し、左辺は流体粒子の密度時間変化率、右辺は流体粒子の体積変化率を表す。すなわち、式（２１）は、流体粒子の密度変化率と体積変化率との関係を表している。
【００５３】
式（２２）は、流体粒子についての運動方程式である。式（２２）の左辺は流体粒子の慣性項、右辺第１項は流体の圧力による力、右辺第２項は重力を表している。
式（２１）〜（２３）をGodunov SPH 法に基づいて離散化し、各粒子の速度に擾乱項を付加すると、以下のような計算式が得られる。
【００５４】
【数１１】
2012111082.tif 000021
【００５５】
式（２４）〜（２７）において、
2012111082.tif 000022
は、それぞれ時刻ｎにおける粒子ａの位置ベクトル、速度ベクトル、及び密度を表す。粒子ｂについても同様の表記が用いられる。ｄｔは時間刻み幅を表し、前の時刻ｎにおける各粒子の位置ベクトル、速度ベクトル、及び密度が既知であるとき、式（２４）〜（２７）により、現在時刻ｎ＋１における粒子ａの位置ベクトル、速度ベクトル、及び密度が求められる。
【００５６】
ここで、Ｗはカーネル関数であり、
2012111082.tif 000023
は再標準化行列である。３次元の場合、再標準化行列は以下のように定義される。
【００５７】
【数１２】
2012111082.tif 000024
【００５８】
式（２４）〜（２８）において、
2012111082.tif 000025
は、以下のように定義される。
【００５９】
【数１３】
2012111082.tif 000026
【００６０】
式（２５）及び（２６）において、
2012111082.tif 000027
は、粒子ａと粒子ｂの間での１次元リーマン問題を解いて得られた時空間の中間値を表し、具体的には以下のように決定される。
まず、粒子ａ及びｂに関して、以下の特性関数を定める。
【００６１】
【数１４】
2012111082.tif 000028
【００６２】
さらに、勾配を以下のように計算する。
【００６３】
【数１５】
2012111082.tif 000029
【００６４】
式（４６）は、速度勾配テンソルを表している。式（４１）〜（５０）を用いて、上記時空間の中間値を以下のように決定する。
【００６５】
【数１６】
2012111082.tif 000030
【００６６】
ところで、式（２５）の右辺第４項が速度に付加された擾乱項であり、以下のように決定される。
時刻ｎにおける粒子ａのエネルギーは、以下のように書くことができる。
【００６７】
【数１７】
2012111082.tif 000031
【００６８】
式（６１）のエネルギーを用いて、擾乱項に含まれる
2012111082.tif 000032
を以下のように定義する。
【００６９】
【数１８】
2012111082.tif 000033
【００７０】
次に、乱数γ1 及びγ2 を用いて、乱数ベクトルを以下のように決定する。
【００７１】
【数１９】
2012111082.tif 000034
【００７２】
乱数γ1 及びγ2 としては、例えば、−１から１までの一様乱数が用いられるが、これ以外の乱数を用いても構わない。そして、粒子ａの速度ベクトルを基に、擾乱の向きを表す単位ベクトルを以下のように決定する。
【００７３】
【数２０】
2012111082.tif 000035
【００７４】
図５は、図４のシミュレーション装置が式（２４）〜（６４）を用いて３次元粒子群の位置を解析するシミュレーション処理の例を示すフローチャートである。
この場合、式（２４）〜（６４）に相当する計算式の情報が格納部４０２に格納されている。計算式の情報には、それぞれの式に含まれる変数、定数等の情報と、それらの変数及び定数を用いた演算の情報が含まれる。処理部４０１は、この計算式の情報を参照しながら、初期時刻０から時刻Ｎ（Ｎは１以上の整数）までの各時刻における各粒子の位置ベクトル、速度ベクトル、及び密度を計算する。
【００７５】
まず、処理部４０１は、時刻を表す変数ｎに０を設定し（ステップ５０１）、初期条件４１１を取得する（ステップ５０２）。初期条件４１１はユーザから与えられる場合もあり、処理部４０１が所定のアルゴリズムで生成する場合もある。この初期条件４１１には、初期時刻０における各粒子の位置ベクトル、速度ベクトル、密度、及びエネルギーの値と、位置ベクトルが取り得る値の範囲等が含まれる。
【００７６】
次に、処理部４０１は、式（６１）により、粒子群中の最初の粒子を処理対象として選択し、時刻ｎにおける処理対象粒子のエネルギーを計算する（ステップ５０３）。そして、時刻ｎにおけるエネルギーを初期時刻０におけるエネルギー（初期エネルギー）と比較する（ステップ５０４）。
【００７７】
時刻ｎにおけるエネルギーが初期エネルギー以下であれば（ステップ５０４，ＹＥＳ）、式（６２）により
2012111082.tif 000036
を計算する（ステップ５０６）。時刻ｎにおけるエネルギーが初期エネルギーより大きければ（ステップ５０４，ＮＯ）、式（６２）に従って
2012111082.tif 000037
を０に設定する（ステップ５０５）。
【００７８】
ステップ５０５の処理を行うことで、処理対象粒子のエネルギーが初期エネルギーを超えた場合に速度の擾乱項の効果をなくすことができ、ランダムな運動エネルギーの過剰な注入が防止される。
【００７９】
次に、処理部４０１は、乱数γ1 及びγ2 を発生し、式（６３）及び（６４）を用いて、それぞれの粒子の速度についての擾乱項（式（２５）の右辺第４項）の値を計算し、その値を格納部４０２に格納する（ステップ５０７）。そして、時刻ｎにおける各粒子の位置ベクトル、速度ベクトル、及び密度と速度の擾乱項の値を用いて、式（２４）〜（５５）により、時刻ｎ＋１における処理対象粒子の位置ベクトル、速度ベクトル、及び密度を計算する（ステップ５０８）。
【００８０】
次に、処理部４０１は、粒子群中のすべての粒子について処理が終了したか否かを判定する（ステップ５０９）。未処理の粒子が残っていれば（ステップ５０９，ＮＯ）、次の粒子を処理対象として選択して、ステップ５０３〜５０８の処理を繰り返す。そして、すべての粒子について処理が終了すれば（ステップ５０９，ＹＥＳ）、時刻ｎ＋１における各粒子の位置ベクトルを出力部４０３に転送する。
【００８１】
出力部４０３は、転送された各粒子の位置ベクトルを、時刻ｎ＋１におけるシミュレーション結果４１２として所定の形式で出力する（ステップ５１０）。シミュレーション結果４１２の出力方法としては、表示装置の画面に粒子群の画像情報として表示したり、有線又は無線の通信回線を介してプリンタや情報処理装置に転送したりする方法が考えられる。シミュレーション結果４１２として、各粒子の速度ベクトル、密度等を含めて出力しても良い。
【００８２】
次に、処理部４０１は、ｎがＮ−１に達したか否かを判定する（ステップ５１１）。ｎがＮ−１より小さければ（ステップ５１１，ＮＯ）、ｎ＝ｎ＋１と設定して（ステップ５１２）、ステップ５０３〜５１０の処理を繰り返す。そして、ｎがＮ−１に達すると（ステップ５１１，ＹＥＳ）、シミュレーション処理を終了する。
【００８３】
なお、図５のシミュレーション処理において、必ずしもすべてのステップを実行する必要はなく、シミュレーション対象や初期条件に応じて一部のステップを省略することも可能である。例えば、ランダムな運動エネルギーが過剰に注入される可能性が低い場合は、ステップ５０４及び５０５の処理を省略しても良い。
【００８４】
３次元粒子群の代わりに２次元粒子群の位置を解析する場合は、３次元ベクトルを２次元ベクトルに変更し、計算式の一部を変更することで、図５と同様のシミュレーション処理を行うことができる。
【００８５】
この場合、式（２８）の再標準化行列の代わりに、以下のような再標準化行列が用いられる。
【００８６】
【数２１】
2012111082.tif 000038
【００８７】
また、式（４６）の速度勾配テンソルの代わりに、以下のような速度勾配テンソルが用いられる。
【００８８】
【数２２】
2012111082.tif 000039
【００８９】
さらに、乱数γ1 及びγ2 の代わりに乱数γを用いて、以下のように回転行列Ａが決定される。
【００９０】
【数２３】
2012111082.tif 000040
【００９１】
乱数γとしては、例えば、−１から１までの一様乱数が用いられるが、これ以外の乱数を用いても構わない。そして、式（６４）の代わりに次式を用いて、速度の擾乱の向きを表す単位ベクトルを決定する。
【００９２】
【数２４】
2012111082.tif 000041
【００９３】
ここで、自由表面流での一般的なベンチマークテストとして、２次元流体のダム崩壊問題に対するシミュレーションの結果を説明する。ダム崩壊問題は、初期時刻０で壁に接して静止している流体の塊が、一定の重力の元で時間の経過とともに崩れていく現象である。
【００９４】
図６は、初期時刻０における流体の２次元形状を示している。図６において、ｘ軸、ｙ軸、及びｘ＝０．５８４（ｍ）の直線６０２がダムの左壁、底、及び右壁にそれぞれ対応し、ｘ＝０〜０．１４６（ｍ）及びｙ＝０〜０．２９２（ｍ）の範囲に長方形６０１の流体が静止している。
【００９５】
時刻０（ｔ＝０）における初期条件としては、長方形６０１内に格子間隔０．００３６５ｍの一様格子を設けて粒子群を配置し、各粒子の初期速度を
2012111082.tif 000042
（ｍ／ｓ）とし、初期密度を
2012111082.tif 000043
とした。密度ρ0 は１０００ｋｇ／ｍ3 、重力加速度ｇは９．８ｍ／ｓ2 、音速ｃは２３．９２ｍ／ｓとした。時間刻み幅ｄｔとしては１．５２６×１０-5ｓを用いた。
【００９６】
図７及び図８は、式（２５）に速度についての擾乱項を付加しない標準的なGodunov SPH 法によるシミュレーション結果の画像を示しており、ｔ＝０〜９．９９４５８（ｓ）の１２個の時刻における流体の形状を表している。
【００９７】
一方、図９及び図１０は、式（２５）の速度についての擾乱項を付加した場合のシミュレーション結果の画像を示しており、ｔ＝０〜９．９９８９８（ｓ）の１２個の時刻における流体の形状を表している。
【００９８】
速度の擾乱項を付加しない場合には、図８のｔ＝６．８７６０２（ｓ）以降において波がほとんど消えてしまっている。これに対して、速度の擾乱項を付加した場合には、図１０のｔ＝６．８７４３２（ｓ）以降においても波が残っており、速度の擾乱項を付加しない場合の図８のｔ＝６．８７６０２（ｓ）以降において見られた波の過剰な減衰が抑止されることにより、波の減衰性が改善されていることが分かる。
【００９９】
上述したように、式（２５）では各粒子の速度に擾乱を与えているが、速度の代わりに位置、加速度、圧力、又は密度に擾乱を与えても良い。例えば、各粒子の位置に擾乱を与える場合は、式（２５）及び（２７）を以下のように書き換えて、図５のシミュレーション処理を実行することができる。
【０１００】
【数２５】
2012111082.tif 000044
【０１０１】
式（７６）は式（２５）から右辺第４項の速度についての擾乱項を取り除いた計算式に相当し、式（７７）の右辺第３項が位置に付加された擾乱項である。この場合、処理部４０１は、図５のステップ５０７において、式（２５）の右辺第４項を計算する代わりに、式（７７）の右辺第３項を計算する。乱数を用いて擾乱の向きを表す単位ベクトルを決定する方法は、速度に擾乱を与える場合と同様である。そして、図５のステップ５０８において、式（２５）及び（２７）の代わりに式（７６）及び（７７）を用いて、時刻ｎ＋１における処理対象粒子の位置ベクトル、速度ベクトル、及び密度を計算する。
【０１０２】
次に、各粒子の加速度に擾乱を与える場合は、式（２５）を以下のように書き換えて、図５のシミュレーション処理を実行することができる。
【０１０３】
【数２６】
2012111082.tif 000045
【０１０４】
式（７８）は式（２５）から右辺第４項の速度についての擾乱項を取り除き、右辺第２項の
2012111082.tif 000046
を
2012111082.tif 000047
に置き換えた計算式に相当する。式（７８）の右辺第２項のうち、
2012111082.tif 000048
を含む部分が加速度に付加された擾乱項である。この場合、処理部４０１は、図５のステップ５０７において、式（２５）の右辺第４項を計算する代わりに、式（７８）の擾乱項を計算する。このとき、乱数γ1 及びγ2 の代わりに乱数γを用いて擾乱項を計算する。そして、図５のステップ５０８において、式（２５）の代わりに式（７８）を用いて、時刻ｎ＋１における処理対象粒子の位置ベクトル、速度ベクトル、及び密度を計算する。
【０１０５】
次に、各粒子の密度に擾乱を与える場合は、式（２５）を式（７６）に置き換え、式（２６）を以下のように書き換えることができる。
【０１０６】
【数２７】
2012111082.tif 000049
【０１０７】
式（７９）の右辺第３項の乱数δが密度に付加された擾乱項である。この場合、処理部４０１は、図５のシミュレーション処理の代わりに、図１１に示すようなシミュレーション処理を実行する。
【０１０８】
図１１のシミュレーション処理は、図５のシミュレーション処理からステップ５０３〜５０６を除いた処理に相当し、ステップ１１０１〜１１０８の処理は、図５のステップ５０１、５０２、及び５０７〜５１２の処理に対応する。
【０１０９】
ただし、処理部４０１は、ステップ１１０３において、式（２５）の右辺第４項を計算する代わりに、乱数δを発生し、その値を密度についての擾乱項の値として格納部４０２に格納する。乱数δとしては、例えば、−０．００００２から０．００００２までの一様乱数が用いられるが、これ以外の乱数を用いても構わない。
【０１１０】
そして、ステップ１１０４において、時刻ｎにおける各粒子の位置ベクトル、速度ベクトル、及び密度と密度の擾乱項の値を用いて、時刻ｎ＋１における処理対象粒子の位置ベクトル、速度ベクトル、及び密度を計算する。この場合、式（２４）、（７６）、（７９）、及び（２７）〜（５５）が計算式として用いられる。
【０１１１】
次に、各粒子の圧力に擾乱を与える場合は、式（２５）を式（７６）に置き換え、圧力の中間値を求める式（５５）を以下のように書き換えて、図１１のシミュレーション処理を実行することができる。
【０１１２】
【数２８】
2012111082.tif 000050
【０１１３】
式（８０）の右辺第２項の乱数δp が圧力に付加された擾乱項である。この場合、処理部４０１は、ステップ１１０３において、乱数δp を発生し、その値を圧力についての擾乱項の値として格納部４０２に格納する。乱数δp としては、例えば、−０．００００２から０．００００２までの一様乱数が用いられるが、これ以外の乱数を用いても構わない。
【０１１４】
そして、ステップ１１０４において、時刻ｎにおける各粒子の位置ベクトル、速度ベクトル、及び密度と圧力の擾乱項の値を用いて、時刻ｎ＋１における処理対象粒子の位置ベクトル、速度ベクトル、及び密度を計算する。この場合、式（２４）、（７６）、（２６）〜（５４）、及び（８０）が計算式として用いられる。
【０１１５】
図４のシミュレーション装置は、水面波の伝播等の波動解析に限らず、流体、弾性体等の連続体や粉粒体の運動を解析する他のシミュレーション処理にも利用できる。また、シミュレーション処理で用いられる計算式は式（２４）〜（８０）に限られるものではなく、各粒子の物理量に擾乱を与える他の計算式を用いても構わない。
【０１１６】
図４のシミュレーション装置は、例えば、図１２に示すような情報処理装置（コンピュータ）を用いて実現することが可能である。図１２の情報処理装置は、Central Processing Unit （ＣＰＵ）１２０１、メモリ１２０２、入力装置１２０３、出力装置１２０４、外部記憶装置１２０５、媒体駆動装置１２０６、及びネットワーク接続装置１２０７を備える。これらはバス１２０８により互いに接続されている。
【０１１７】
メモリ１２０２は、例えば、Read Only Memory（ＲＯＭ）、Random Access Memory（ＲＡＭ）、フラッシュメモリ等の半導体メモリであり、シミュレーション処理に用いられるプログラム及びデータを格納する。例えば、ＣＰＵ１２０１は、メモリ１２０２を利用してプログラムを実行することにより、図５又は図１１のシミュレーション処理を行う。メモリ１２０２は、図４の格納部４０２としても使用できる。
【０１１８】
入力装置１２０３は、例えば、キーボード、ポインティングデバイス等であり、オペレータからの指示や情報の入力に用いられる。出力装置１２０４は、例えば、表示装置、プリンタ、スピーカ等であり、オペレータへの問い合わせや処理結果の出力に用いられる。出力装置１２０４は、図４の出力部４０３としても使用できる。
【０１１９】
外部記憶装置１２０５は、例えば、磁気ディスク装置、光ディスク装置、光磁気ディスク装置、テープ装置等である。この外部記憶装置１２０５には、ハードディスクドライブも含まれる。情報処理装置は、この外部記憶装置１２０５にプログラム及びデータを格納しておき、それらをメモリ１２０２にロードして使用することができる。
【０１２０】
媒体駆動装置１２０６は、可搬型記録媒体１２０９を駆動し、その記録内容にアクセスする。可搬型記録媒体１２０９は、メモリデバイス、フレキシブルディスク、光ディスク、光磁気ディスク等である。この可搬型記録媒体１２０９には、Compact Disk Read Only Memory （ＣＤ−ＲＯＭ）、Digital Versatile Disk（ＤＶＤ）、Universal Serial Bus（ＵＳＢ）メモリ等も含まれる。オペレータは、この可搬型記録媒体１２０９にプログラム及びデータを格納しておき、それらをメモリ１２０２にロードして使用することができる。
【０１２１】
このように、シミュレーション処理に用いられるプログラム及びデータを格納するコンピュータ読み取り可能な記録媒体には、メモリ１２０２、外部記憶装置１２０５、及び可搬型記録媒体１２０９のような、物理的な（非一時的な）記録媒体が含まれる。
【０１２２】
ネットワーク接続装置１２０７は、通信ネットワーク１０３に接続され、通信に伴うデータ変換を行う通信インタフェースである。情報処理装置は、プログラム及びデータを外部の装置からネットワーク接続装置１２０７を介して受け取り、それらをメモリ１２０２にロードして使用することができる。ネットワーク接続装置１２０７は、図４の出力部４０３としても使用できる。
【０１２３】
開示の実施形態とその利点について詳しく説明したが、当業者は、特許請求の範囲に明確に記載した本発明の範囲から逸脱することなく、様々な変更、追加、省略をすることができるであろう。
【図１】
2012111082.tif 000051
【図２】
2012111082.tif 000052
【図３】
2012111082.tif 000053
【図４】
2012111082.tif 000054
【図５】
2012111082.tif 000055
【図６】
2012111082.tif 000056
【図７】
2012111082.tif 000057
【図８】
2012111082.tif 000058
【図９】
2012111082.tif 000059
【図１０】
2012111082.tif 000060
【図１１】
2012111082.tif 000061
【図１２】
2012111082.tif 000062
【手続補正書】
【提出日】平成25年8月22日(2013.8.22)
【手続補正１】
【補正対象書類名】特許請求の範囲
【補正対象項目名】全文
【補正方法】変更
【補正の内容】
【特許請求の範囲】
【請求項１】
格納部に格納された、現在時刻における複数の粒子の位置、速度、加速度、圧力、又は密度に付加された擾乱項の情報を用いて、該擾乱項の値を求め、
前記複数の粒子の前記擾乱項の値を前記格納部に格納し、
前記格納部に格納された、前記複数の粒子の前記擾乱項を含む計算式であって、前の時刻における該複数の粒子の位置、速度、及び密度から前記現在時刻における該複数の粒子の位置、速度、及び密度を求める該計算式の情報を用い、初期時刻における該複数の粒子の位置、速度、及び密度を該前の時刻における該複数の粒子の位置、速度、及び密度として用い、前記複数の粒子の前記擾乱項の値を該計算式に適用することで、該現在時刻における該複数の粒子の位置、速度、及び密度を求め、
得られた前記複数の粒子の位置、速度、及び密度を前記前の時刻における前記複数の粒子の位置、速度、及び密度として用いて該複数の粒子の前記擾乱項の値と前記現在時刻における該複数の粒子の位置、速度、及び密度を求める処理を繰り返すことで、前記初期時刻より後の複数の時刻における該複数の粒子の位置、速度、及び密度を順次求め、
前記複数の時刻における前記複数の粒子の位置をシミュレーション結果として出力する
処理をコンピュータに実行させるプログラム。
【請求項２】
前記コンピュータは、前記現在時刻における各粒子の位置又は速度に付加された擾乱の向きを示す乱数を用いて、前記複数の粒子の位置又は速度に付加された前記擾乱項の値を求めることを特徴とする請求項１記載のプログラム。
【請求項３】
前記コンピュータは、前記現在時刻における各粒子の加速度、圧力、又は密度に付加された擾乱の大きさを示す乱数を用いて、前記複数の粒子の加速度、圧力、又は密度に付加された前記擾乱項の値を求めることを特徴とする請求項１記載のプログラム。
【請求項４】
前記コンピュータは、前記前の時刻における各粒子の位置、速度、及び密度を用いて該前の時刻における各粒子のエネルギーを求め、得られたエネルギーが前記初期時刻における該粒子のエネルギーより大きい場合、該粒子の擾乱項の値を零に設定することを特徴とする請求項１、２、又は３記載のプログラム。
【請求項５】
現在時刻における複数の粒子の位置、速度、加速度、圧力、又は密度に付加された擾乱項を含み、前の時刻における該複数の粒子の位置、速度、及び密度から該現在時刻における該複数の粒子の位置、速度、及び密度を求める計算式の情報を格納する格納部と、
前記複数の粒子の前記擾乱項の値を求め、該複数の粒子の該擾乱項の値を前記格納部に格納し、初期時刻における該複数の粒子の位置、速度、及び密度を前記前の時刻における前記複数の粒子の位置、速度、及び密度として用い、該複数の粒子の該擾乱項の値を前記計算式に適用することで、前記現在時刻における該複数の粒子の位置、速度、及び密度を求め、得られた前記複数の粒子の位置、速度、及び密度を前記前の時刻における前記複数の粒子の位置、速度、及び密度として用いて該複数の粒子の該擾乱項の値と該現在時刻における該複数の粒子の位置、速度、及び密度を求める処理を繰り返すことで、該初期時刻より後の複数の時刻における該複数の粒子の位置、速度、及び密度を順次求める処理部と、
前記複数の時刻における前記複数の粒子の位置をシミュレーション結果として出力する出力部と
を備えることを特徴とするシミュレーション装置。
【請求項６】
コンピュータにより実行されるシミュレーション方法であって、
格納部に格納された、現在時刻における複数の粒子の位置、速度、加速度、圧力、又は密度に付加された擾乱項の情報を用いて、該擾乱項の値を求め、
前記複数の粒子の前記擾乱項の値を前記格納部に格納し、
前記格納部に格納された、前記複数の粒子の前記擾乱項を含む計算式であって、前の時刻における該複数の粒子の位置、速度、及び密度から前記現在時刻における該複数の粒子の位置、速度、及び密度を求める該計算式の情報を用い、初期時刻における該複数の粒子の位置、速度、及び密度を該前の時刻における該複数の粒子の位置、速度、及び密度として用い、前記複数の粒子の前記擾乱項の値を該計算式に適用することで、該現在時刻における該複数の粒子の位置、速度、及び密度を求め、
得られた前記複数の粒子の位置、速度、及び密度を前記前の時刻における前記複数の粒子の位置、速度、及び密度として用いて該複数の粒子の前記擾乱項の値と前記現在時刻における該複数の粒子の位置、速度、及び密度を求める処理を繰り返すことで、前記初期時刻より後の複数の時刻における該複数の粒子の位置、速度、及び密度を順次求め、
前記複数の時刻における前記複数の粒子の位置をシミュレーション結果として出力する
ことを特徴とするシミュレーション方法。
【国際調査報告】
2012111082.tif 000063
2012111082.tif 000064
2012111082.tif 000065
2012111082.tif 000066