JP6163897B2 - 数値計算プログラム、数値計算方法及び情報処理装置 - Google Patents

Description

本発明は、数値計算プログラム、数値計算方法及び情報処理装置に関する。

流体や弾性体等の連続体の振る舞いを解く数値計算手法としては、格子をベースにして微分方程式の近似解を求める有限差分法、有限要素法、有限体積法等が知られている。近年では数値計算をＣＡＥ（Computer Aided Engineering）などの分野で活用するため、これらの数値計算手法も発展し、流体と構造物が相互作用する問題が解かれるようになってきた。しかしながら、格子を用いる手法では、自由表面などの界面の存在する問題や、流体−構造連成問題などの移動境界が発生する場合には取り扱いが複雑なため、プログラム作成が困難である場合が多い。

これに対して、格子を用いない粒子法（ＭＰＳ（Moving Particle Semi-implicit）法、ＳＰＨ（Smoothed Particle Hydrodynamics）法等）では、移動境界の取り扱いに特別な処置を必要としない。そのため近年、粒子法は広く用いられるようになってきている。

また、粒子法により高精度に計算する手法として、リーマンソルバを用いた計算手法が知られている。ある文献には、粒子法の計算で問題となる圧力の振動を抑制することで、精度よくシミュレーションを行う技術が開示されている。

例えば、１次元の非粘性流体の方程式におけるリーマン問題解を利用する計算手法（リーマンソルバと呼ぶ）を導入することにより、計算の高精度化を行うことを開示した文献が存在する。

ここで、リーマンソルバを粒子法に導入する方法の一例について簡単に説明する。僅かに圧縮する流体の方程式（０−１）乃至（０−３）を考える。

式（０−１）は質量保存則、式（０−２）が運動量保存則、式（０−３）は状態方程式である。ここでｖ（太字ｖ）は流体の速度ベクトル、ρは密度、ｐは圧力、ρ₀は圧力が０となる密度（すなわち基準密度）である。また、ｃは音速であり、ｔは時刻を表す。

そして、式（０−１）乃至（０−３）を以下のように離散化する。

ここで、ｘ（太字ｘ）は３次元空間の位置ベクトルを表す。また、ｘ_i、ｖ_i、ρ_i及びｍ_iは、それぞれ、粒子ｉの位置ベクトル、速度ベクトル、密度、質量を表す。ｘ_ij ^n+1/2は、粒子ｉと粒子ｊの相対位置ベクトルである。具体的には、以下のように表される。
ｘ_ij ^n+1/2＝ｘ_i ^n+1/2−ｘ_j ^n+1/2

また、上付き添え字はシミュレーションステップでの時刻を表わす。また、式（０−４）乃至（０−７）に含まれる記号は以下のように定義される。

さらに、Ｗはカーネル関数（重み関数とも呼ぶ）であり、以下のスプライン関数が良く用いられる。

ここでｈは、粒子間の影響半径で、初期状態の平均粒子間隔の２倍から３倍程度が良く用いられる。βはカーネル関数の全空間積分量が１になるように調整された値で、２次元の場合は０．７πｈ² 、３次元の場合はπｈ³と決められる。

式（０−８）及び（０−９）は、粒子ｉｊ間の中間値であって、１次元の流体の方程式を数値的に解き、ｄｔ／２だけ進んだ時刻における解を用いる。すなわち、図１に示すように、粒子ｉｊ間に軸ｌ_ijをセットして、この軸上のベクトルに射影した速度ベクトルは、それぞれ以下のように表される。

具体的な決定方法は以下のようなものである。まず、式（０−１）乃至（０−３）について、空間次元が１次元の場合の方程式を整理する。

そして、式（０−１１）についてリーマン不変量ｑ⁺及びｑ^-を用いて特性方程式の形式に書き直す。

そして、時刻ｎにおける粒子ｉ及びｊについてｑ⁺及びｑ^-を以下のように表す。

さらに、空間勾配を以下のように求めておく。

このようにして変数を整えた後に、図２に示すように、粒子ｉｊ間に物理量が不連続な面が存在する状況を初期値とした問題（リーマン問題）を考え、式（０−１２）及び（０−１３）の解析解を用いてｑ⁺及びｑ^-の粒子ｉｊ間の中間における、時刻ｎ＋１／２の値を予測する。具体的には以下のように算出される。

このような式（０−１６）及び（０−１７）を、式（０−５）及び（０−６）に代入することで粒子法のシミュレーションを高精度に実行できる。

一方、剛体と流体の連成計算や、混ざることのない別々の液体（たとえば水と油）の混流現象のシミュレーションは、自由表面や移動境界の扱いが簡単な粒子法が期待されている。混流現象の場合は安定状態となるときの密度が違う液体を取り扱うため、標準的なＳＰＨ法だけでは一般に計算は難しくなる。

密度の違う液体の混流を粒子法で扱う計算方法として、ある文献では、以下のような計算を行うことを示している。すなわち、密度の異なる流体方程式として以下の式（０−１８）乃至（０−２１）を用いる。

ｇは重力による加速度であり、ρ_sは基準密度である。式（０−１）乃至（０−３）と異なる点は基準密度が位置によって異なることである。基準密度の空間変化で、異なる流体の混流を表している。

そして、式（０−１８）乃至（０−２１）をＳＰＨ法で以下のように各粒子について離散化する。

式（０−２２）は質量保存則、式（０−２３）は運動量保存則、式（０−２４）は状態方程式を表す。ζ及びηはパラメータで定数である。ここでρ_s,iは粒子ｉの基準密度を表す。
ｘ_ij＝ｘ_i−ｘ_j
ｖ_ij＝ｖ_i−ｖ_j

式（０−２２）乃至（０−２４）について、一般的な常微分方程式の数値解析法であるオイラー法やリープフロッグ法などで時間発展を計算することでシミュレーションが可能となる。

しかしながら、粒子法で用いられているリーマンソルバによる計算の高精度化は、この技術の基本的な式（０−１８）乃至（０−２１）には適用することができない。従って、密度差がある流体の混流に対して計算精度が低下する。

国際公開第２０１２／１１１０８２号パンフレット

Shu-ichiro Inutsuka, "Reformulation of Smoothed Particle Hydrodynamics with Riemann Solver", Journal of Computaional Physics, 179, (2002) 238-267 Paul W. Cleary, "Extension of SPH to predict feeding freezing and defect creation in low pressure die casting", Applied Mathematical Modelling, 34(2010), pp. 3189-3201.

従って、本発明の目的は、基準密度が異なる流体の混流現象を精度良く数値計算する技術を提供することである。

本発明に係る数値計算方法は、（Ａ）（ａ１）第１の基準密度の流体に係る第１の粒子の密度の代わりに、当該第１の粒子の密度の、第１の基準密度に対する比率をリーマン不変量において用いることで得られる第１の物理量を算出し、（ａ２）第２の基準密度の流体に係る第２の粒子の密度の代わりに、当該第２の粒子の密度の、第２の基準密度に対する比率をリーマン不変量において用いることで得られる第２の物理量を算出し、（ａ３）第１の物理量の勾配及び第２の物理量の勾配を算出し、（ａ４）第１の粒子と第２の粒子との間における第１及び第２の物理量に関する時空間中間値を、第１の物理量、第２の物理量、第１の物理量の勾配及び第２の物理量の勾配を用いて算出し、（ａ５）第１の粒子と第２の粒子との間における第１及び第２の物理量に関する中間値から、第１の粒子と第２の粒子との間における圧力の時空間中間値を算出し、（ａ６）第１の粒子と第２の粒子との間における第１及び第２の物理量に関する中間値から、第１の粒子と第２の粒子との間における速度の時空間中間値を算出する第１の処理を、第１の粒子との位置関係から特定される第２の粒子の各々について実行し、（Ｂ）算出された圧力の時空間中間値を用いて第１の粒子の速度を更新し、（Ｃ）算出された速度の時空間中間値に基づき第１の粒子の密度を更新する処理を含む。

一側面によれば、基準密度が異なる流体の混流現象を精度良く数値計算できるようになる。

図１は、粒子ｉｊの中間値を説明するための図である。 図２は、粒子ｉｊ間に物理量が不連続な面が存在する状況を表す図である。 図３は、本実施の形態において想定される混流の例を示す図である。 図４は、本実施の形態に係る情報処理装置の機能ブロック図である。 図５は、本実施の形態に係る処理の処理フローを示す図である。 図６は、本実施の形態に係る処理の処理フローを示す図である。 図７は、本実施の形態に係る処理の処理フローを示す図である。 図８は、本実施の形態に係る処理の処理フローを示す図である。 図９は、本実施の形態に係る処理の処理フローを示す図である。 図１０は、本実施の形態に係る処理の処理フローを示す図である。 図１１は、コンピュータの機能ブロック図である。

まず、本発明の実施の形態に係る演算の考え方について説明する。

本実施の形態では、以下の流体運動を表す方程式系を粒子法で解く。

ここでｖ及びｐは流体の速度場及び圧力場を表す。また、ρ及びρ_sは流体の密度及び基準密度（圧力が０の状態となる密度）を表す。また、ｃは音速を表す。

式（１）乃至（４）は、図３に示すような、水と油のように互いが混じり合うことなく、密度が異なる流体の混流を表す場合などに用いられる。

この式（１）乃至（４）を粒子法で空間離散化すると、以下のように表される。

なお、ｊは例えば影響範囲内の他の粒子を表す。

さらに、時間について離散化することで以下の漸化式を導出する。

ここで、ｄｔは時間刻み、その他の記号は以下のような意味である。

式（１０）のように、密度の時間発展の右辺第二項におけるρⁿ _i／ρⁿ _jの部分以外は、時刻ｎ＋１／２の値と時刻ｎとｎ＋１の平均値とを含むようになる。そうするとρⁿ _i／ρⁿ _jの時間変化が小さい状況では、２次の精度を持つような計算方法となる。

ｖ_ij ^lij,n+1/2及びｐ_ij ^lij,n+1/2は粒子ｉｊ間の中間値であり、以下の手順で算出する。

まず、式（１）乃至（４）について１次元の方程式を考える。

ここで、リーマン不変量におけるρに代わってρ／ρ_sを用いてリーマン不変量を変形することで、以下のような物理量を新たに導入する。

そして、これらの物理量を用いて整理すると、以下のような式が得られる。

式（１１−１）及び（１１−２）の右辺第一項だけを考えると、式（０−１２）及び（０−１３）と同様の形式であるため、右辺の第一項のみを考えて時間発展を以下のように行い、右辺第二項及び第三項を補う項を後で追加するようにして計算を行う。
Ａ）より具体的には、変数ｑ⁺及びｑ^-を以下のようにして求める。

これらは、時刻ｎにおける粒子ｉ及びｊについての上記物理量である。
Ｂ）また、変数ｑ⁺及びｑ^-の勾配を以下のように計算する。

式（１６）乃至（１９）における記号は以下のように表される。

（２１）式は、変数ｑ⁺及びｑ^-の勾配が、他の粒子の基準密度との関係に基づく基準密度の空間勾配を用いて算出されることを表している。すなわち、この分が、背景技術で述べた混流ではない場合との差である。
Ｃ）変数ｑ⁺及びｑ^-について、粒子ｉｊ間の中間値を以下のように計算する。

式（２３）及び（２４）の右辺第二項はそれぞれ式（０−１４）及び（０−１５）の右辺第二項と類似の形式であり、リーマンソルバの考え方を用いて得られた項であることが分かる。

なお、ρ_ij ⁿ及びρ_s,ij ⁿは粒子ｉｊ間での密度、及び基準密度の中間値で、以下のように表される。

式（２３）及び（２４）の右辺第三項及び第四項は、それぞれ式（１１−１）及び（１１−２）の右辺第二項及び第三項を差分近似したものである。式（２３）式（２４）の右辺第三項のDρ_s /Dt|_i およびDρ_s /Dt|_j は、０として問題ない。また、基準密度の空間微分は、以下のように表され、標準的なＳＰＨ法の計算手法により得られる。

Ｄ）上で算出した中間値ｑ_ij ^n+1/2,+及びｑ_ij ^n+1/2,-を用いて式（２７）及び（２８）より速度及び圧力の中間値を算出して、式（９）及び（１０）に代入する。

以上の考え方に基づき、具体的な処理及び当該処理を行うための装置について説明する。

図４に、本実施の形態に係る情報処理装置１００の機能ブロック図を示す。情報処理装置１００は、初期条件データ格納部１１０と、処理部１２０と、データ格納部１３０と、処理結果格納部１４０と、出力部１５０とを有する。そして、情報処理装置１００は、出力装置２００に接続されている。

初期条件データ格納部１１０には、初期位置、初期速度、初期密度、基準密度、時間刻みｄｔ、粒子にかかる外力（例えば重力）のデータなどが格納されている。

処理部１２０は、初期条件データ格納部１１０に格納されているデータを用いて、各タイムステップについて、各粒子についての速度、密度及び位置等のデータを算出し、処理結果格納部１４０に格納する。なお、処理途中のデータについては、データ格納部１３０に格納される。

出力部１５０は、処理結果格納部１４０に格納されるデータを、出力装置２００（例えば印刷装置、表示装置、場合によってはネットワークで接続される他のコンピュータなどの装置）に出力する。

次に、図５乃至図１０を用いて、情報処理装置１００の処理内容について説明する。

まず、処理部１２０は、初期条件データ格納部１１０に格納されている初期条件データを読み出す（図５：ステップＳ１）。また、処理部１２０は、ｎ＝１と設定する（ステップＳ３）。

そして、処理部１２０は、初期条件データにおいて設定がなされている粒子のうち未処理の粒子ｉを１つ特定する（ステップＳ５）。

その後、処理部１２０は、粒子ｉの位置をｄｔ／２だけ更新する（ステップＳ７）。具体的には、以下のような演算を実行する。ｘ_i ¹については、初期条件データに含まれる。

また、処理部１２０は、粒子ｉの加速度、密度の時間変化項を初期化する（ステップＳ９）。具体的には、以下のように設定を行う。

さらに、処理部１２０は、粒子ｉにかかる外力を算出する（ステップＳ１１）。具体的には、以下のような演算を実行する。外力ｆⁿ _iは、初期条件データに含まれる。

また、処理部１２０は、粒子ｉの影響範囲に入る他の粒子ｊを抽出する（ステップＳ１３）。そして処理は端子Ａを介して図６の処理に移行する。

端子Ａを介して図６の処理の説明に移行して、処理部１２０は、抽出された粒子ｊの中から未処理の粒子ｊを１つ特定する（ステップＳ１５）。そして、処理部１２０は、粒子ｉｊ間の圧力に関する時空間中間値の算出処理を実行する（ステップＳ１７）。この処理については、図７及び図８を用いて説明する。

処理部１２０は、物理量ｑの中間値算出処理を実行する（図７：ステップＳ３１）。具体的には、図８の処理を実行する。

処理部１２０は、粒子ｉｊの密度、基準密度及び速度から、式（１３）及び（１４）に従って、物理量ｑ^- _i、ｑ⁺ _jを算出する（ステップＳ３５）。

また、処理部１２０は、粒子ｉｊの密度、基準密度及び速度の勾配から、式（１７）及び（１８）に従って、物理量ｑの勾配を算出する（ステップＳ３７）。

さらに、処理部１２０は、物理量ｑの時空間中間値ｑ_ij ^n+1/2,+及びｑ_ij ^n+1/2,-を、式（２３）及び（２４）に従って算出する（ステップＳ３９）。

図８の演算結果については、データ格納部１３０の容量が十分にあれば保持しておくことで、後の処理を省略できる。データ格納部１３０の容量が十分なければ、後の処理でも同じ演算を実行する。

図７の処理の説明に戻って、処理部１２０は、物理量ｑの時空間中間値を用いて、式（２８）に従って、圧力に関する時空間中間値ｐ_ij ^n+1/2を算出する（ステップＳ３３）。

圧力に関する時空間中間値ｐ_ij ^n+1/2を速度に関する時空間中間値とは別に先行して算出するのは、速度ｖ_i ⁿ⁺¹が、速度に関する時空間中間値の計算に用いられるためである。

図６の処理の説明に戻って、処理部１２０は、粒子ｉの加速度を、圧力に関する時空間中間値ｐ_ij ^n+1/2を用いて更新する（ステップＳ１９）。このステップにおける演算は、以下のように表される。

その後、処理部１２０は、ステップＳ１３で抽出された粒子のうち未処理の粒子ｊが存在しているか判断する（ステップＳ２１）。未処理の粒子ｊが存在する場合には、処理はステップＳ１５へ戻る。一方、未処理の粒子ｊが存在しない場合には、処理部１２０は、粒子ｉの加速度から粒子ｉの速度を更新する（ステップＳ２３）。このステップにおける演算は、以下のように表される。

そして、処理部１２０は、ステップＳ１３で抽出された粒子ｊを未処理に設定する（ステップＳ２５）。なお、ステップＳ１３と同様の処理を実行するようにしても良い。処理は、端子Ｂを介して図９の処理に移行する。

図９の処理の説明に移行して、処理部１２０は、ステップＳ１３で抽出された粒子のうち未処理の粒子ｊを１つ特定する（ステップＳ４１）。そして、処理部１２０は、粒子ｉｊ間の速度に関する時空間中間値の算出処理を実行する（ステップＳ４３）。この処理については、図１０を用いて説明する。

処理部１２０は、物理量ｑの中間値算出処理を実行する（ステップＳ６３）。この処理は、図８の処理と同じである。上でも述べたように、物理量ｑの中間値については、データ格納部１３０に格納されている場合には、この処理を実行する代わりに、データ格納部１３０から読み出しても良い。

そして、処理部１２０は、物理量ｑの時空間中間値を用いて、式（２７）に従って、速度に関する時空間中間値ｖ_ij ^lij,n+1/2を算出する（ステップＳ６５）。なお、式（２７）の後に以下の演算を行う。

図９の処理の説明に戻って、処理部１２０は、粒子ｉの密度時間変化項を、以下の式に従って更新する（ステップＳ４５）。

そして、処理部１２０は、ステップＳ１３で抽出された粒子のうち未処理の粒子ｊが存在するか判断する（ステップＳ４７）。未処理の粒子ｊが存在する場合には、処理はステップＳ４１に戻る。一方、未処理の粒子ｊが存在しない場合には、処理部１２０は、粒子ｉの密度時間変化項から粒子ｉの密度を、以下の式に従って更新する（ステップＳ４９）。

さらに、処理部１２０は、粒子ｉの位置を、ステップＳ２３で算出された時刻ｎ＋１における速度ｖⁿ⁺¹ _iを用いてｄｔ／２だけ更新する（ステップＳ５１）。以下に示す式で算出する。

そして、処理部１２０は、未処理の粒子ｉが存在するか判断する（ステップＳ５３）。未処理の粒子ｉが存在する場合には、端子Ｃを介して図５のステップＳ５に戻る。一方、未処理の粒子ｉが存在しない場合には、処理部１２０は、時刻ｎ＋１の結果を、処理結果格納部１４０に格納し、出力部１５０は、処理結果格納部１４０に格納された時刻ｎ＋１の処理結果を出力装置２００に対して所定の形式で出力する（ステップＳ５５）。

そして、処理部１２０は、時刻ｎ＋１は終了時刻であるか判断する（ステップＳ５７）。終了時刻でなければ、処理部１２０は、ｎ＝ｎ＋１と設定し（ステップＳ５９）、粒子ｉを全て未処理に設定する（ステップＳ６１）。そして、処理は端子Ｃを介して図５のステップＳ５に戻る。一方、時刻ｎ＋１が終了時刻であれば、処理を終了する。

以上のような処理を実行することで、粒子法で用いられているリーマンソルバによる計算の高精度化が実現する。

以上本発明の実施の形態を説明したが、本発明はこれに限定されるものではない。例えば、図４の機能ブロック図は一例であって、プログラムモジュール構成及びファイル構成とは一致しない場合もある。

処理フローについても、処理結果が変わらない限り、ステップの処理順番を入れ替えたり、複数ステップを並列に実行するようにしても良い。

さらに、情報処理装置１００は、１台のコンピュータではなく、複数のコンピュータで構成される場合もある。また、クライアントサーバ型システムで実装しても良いし、スタンドアロン型システムで実装しても良い。

なお、上で述べた情報処理装置１００は、コンピュータ装置であって、図１１に示すように、メモリ２５０１とＣＰＵ（Central Processing Unit）２５０３とハードディスク・ドライブ（ＨＤＤ：Hard Disk Drive）２５０５と表示装置２５０９に接続される表示制御部２５０７とリムーバブル・ディスク２５１１用のドライブ装置２５１３と入力装置２５１５とネットワークに接続するための通信制御部２５１７とがバス２５１９で接続されている。オペレーティング・システム（ＯＳ：Operating System）及び本実施例における処理を実施するためのアプリケーション・プログラムは、ＨＤＤ２５０５に格納されており、ＣＰＵ２５０３により実行される際にはＨＤＤ２５０５からメモリ２５０１に読み出される。ＣＰＵ２５０３は、アプリケーション・プログラムの処理内容に応じて表示制御部２５０７、通信制御部２５１７、ドライブ装置２５１３を制御して、所定の動作を行わせる。また、処理途中のデータについては、主としてメモリ２５０１に格納されるが、ＨＤＤ２５０５に格納されるようにしてもよい。本技術の実施例では、上で述べた処理を実施するためのアプリケーション・プログラムはコンピュータ読み取り可能なリムーバブル・ディスク２５１１に格納されて頒布され、ドライブ装置２５１３からＨＤＤ２５０５にインストールされる。インターネットなどのネットワーク及び通信制御部２５１７を経由して、ＨＤＤ２５０５にインストールされる場合もある。このようなコンピュータ装置は、上で述べたＣＰＵ２５０３、メモリ２５０１などのハードウエアとＯＳ及びアプリケーション・プログラムなどのプログラムとが有機的に協働することにより、上で述べたような各種機能を実現する。

以上述べた本実施の形態をまとめると、以下のようになる。

本実施の形態に係る数値計算方法は、（Ａ）（ａ１）第１の基準密度の流体に係る第１の粒子の密度の代わりに、当該第１の粒子の密度の、第１の基準密度に対する比率をリーマン不変量において用いることで得られる第１の物理量を算出し、（ａ２）第２の基準密度の流体に係る第２の粒子の密度の代わりに、当該第２の粒子の密度の、第２の基準密度に対する比率をリーマン不変量において用いることで得られる第２の物理量を算出し、（ａ３）第１の物理量の勾配及び第２の物理量の勾配を算出し、（ａ４）第１の粒子と第２の粒子との間における第１及び第２の物理量に関する時空間中間値を、第１の物理量、第２の物理量、第１の物理量の勾配及び第２の物理量の勾配を用いて算出し、（ａ５）第１の粒子と第２の粒子との間における第１及び第２の物理量に関する中間値から、第１の粒子と第２の粒子との間における圧力の時空間中間値を算出し、（ａ６）第１の粒子と第２の粒子との間における第１及び第２の物理量に関する中間値から、第１の粒子と第２の粒子との間における速度の時空間中間値を算出する第１の処理を、第１の粒子との位置関係から特定される第２の粒子の各々について実行し、（Ｂ）算出された圧力の時空間中間値を用いて第１の粒子の速度を更新し、（Ｃ）算出された速度の時空間中間値に基づき第１の粒子の密度を更新する処理を含む。

このようにリーマン不変量を変形して新たな物理量を定義することで、リーマンソルバによる高精度計算が実現される。

また、上で述べた勾配を算出する処理が、（ａ３１）第１の基準密度の空間勾配を用いて、第１の物理量の勾配を算出し、（ａ３２）第２の基準密度の空間勾配を用いて、第２の物理量の勾配を算出する処理を含むようにしても良い。リーマン不変量を、基準密度に対する粒子の密度の比率を用いて変形することで、リーマン不変量の勾配についても、第１の基準密度の空間勾配を用いるような形に変形される。

さらに、上で述べた第１の粒子の速度を更新する処理が、（ｂ１）圧力の時空間中間値から、第１の粒子の加速度を更新し、（ｂ２）第１の粒子の加速度を用いて、第１の粒子の速度を更新する処理を含むようにしても良い。加速度は、重力などの外力によって初期的に設定されるが、第２の粒子からの影響を重ね合わせて、さらに更新される。

また、上で述べた第１の粒子の密度を更新する処理が、（ｃ１）第１の粒子の速度と速度の時空間中間値とを用いて、第１の粒子の密度時間変化を算出し、（ｃ２）第１の粒子の密度時間変化を用いて、第１の粒子の密度を更新する処理を含むようにしても良い。

例えば、第１の処理の中で第１の粒子との位置関係から特定された第２の粒子について圧力に関する時空間中間値が得られて第１の粒子の速度が更新された後に、第１の粒子の密度時間変化を算出するようにしても良い。このようにすれば計算の精度が向上する。

なお、上で述べたような処理をコンピュータに実行させるためのプログラムを作成することができ、当該プログラムは、例えばフレキシブル・ディスク、ＣＤ−ＲＯＭなどの光ディスク、光磁気ディスク、半導体メモリ（例えばＲＯＭ）、ハードディスク等のコンピュータ読み取り可能な記憶媒体又は記憶装置に格納される。なお、処理途中のデータについては、ＲＡＭ等の記憶装置に一時保管される。

以上の実施例を含む実施形態に関し、さらに以下の付記を開示する。

（付記１）
第１の基準密度の流体に係る第１の粒子の密度の代わりに、当該第１の粒子の密度の、前記第１の基準密度に対する比率をリーマン不変量において用いることで得られる第１の物理量を算出し、
第２の基準密度の流体に係る第２の粒子の密度の代わりに、当該第２の粒子の密度の、前記第２の基準密度に対する比率をリーマン不変量において用いることで得られる第２の物理量を算出し、
前記第１の物理量の勾配及び前記第２の物理量の勾配を算出し、
前記第１の粒子と前記第２の粒子との間における前記第１及び第２の物理量に関する時空間中間値を、前記第１の物理量、前記第２の物理量、前記第１の物理量の勾配及び前記第２の物理量の勾配を用いて算出し、
前記第１の粒子と前記第２の粒子との間における前記第１及び第２の物理量に関する中間値から、前記第１の粒子と前記第２の粒子との間における圧力の時空間中間値を算出し、
前記第１の粒子と前記第２の粒子との間における前記第１及び第２の物理量に関する中間値から、前記第１の粒子と前記第２の粒子との間における速度の時空間中間値を算出する
第１の処理を、前記第１の粒子との位置関係から特定される第２の粒子の各々について実行し、
算出された前記圧力の時空間中間値を用いて前記第１の粒子の速度を更新し、
算出された前記速度の時空間中間値に基づき前記第１の粒子の密度を更新する
処理を、コンピュータに実行させるための数値計算プログラム。

（付記２）
前記勾配を算出する処理が、
前記第１の基準密度の空間勾配を用いて、前記第１の物理量の勾配を算出し、
前記第２の基準密度の空間勾配を用いて、前記第２の物理量の勾配を算出する
処理を含む付記１記載の数値計算プログラム。

（付記３）
前記第１の粒子の速度を更新する処理が、
前記圧力の時空間中間値から、前記第１の粒子の加速度を更新し、
前記第１の粒子の加速度を用いて、前記第１の粒子の速度を更新する
処理を含む付記１又は２記載の数値計算プログラム。

（付記４）
前記第１の粒子の密度を算出する処理が、
前記第１の粒子の速度と前記速度の時空間中間値とを用いて、前記第１の粒子の密度時間変化を算出し、
前記第１の粒子の密度時間変化を用いて、前記第１の粒子の密度を更新する
処理を含む付記１乃至３のいずれか１つ記載の数値計算プログラム。

（付記５）
第１の基準密度の流体に係る第１の粒子の密度の代わりに、当該第１の粒子の密度の、前記第１の基準密度に対する比率をリーマン不変量において用いることで得られる第１の物理量を算出し、
第２の基準密度の流体に係る第２の粒子の密度の代わりに、当該第２の粒子の密度の、前記第２の基準密度に対する比率をリーマン不変量において用いることで得られる第２の物理量を算出し、
前記第１の物理量の勾配及び前記第２の物理量の勾配を算出し、
前記第１の粒子と前記第２の粒子との間における前記第１及び第２の物理量に関する時空間中間値を、前記第１の物理量、前記第２の物理量、前記第１の物理量の勾配及び前記第２の物理量の勾配を用いて算出し、
前記第１の粒子と前記第２の粒子との間における前記第１及び第２の物理量に関する中間値から、前記第１の粒子と前記第２の粒子との間における圧力の時空間中間値を算出し、
前記第１の粒子と前記第２の粒子との間における前記第１及び第２の物理量に関する中間値から、前記第１の粒子と前記第２の粒子との間における速度の時空間中間値を算出する
第１の処理を、前記第１の粒子との位置関係から特定される第２の粒子の各々について実行し、
算出された前記圧力の時空間中間値を用いて前記第１の粒子の速度を更新し、
算出された前記速度の時空間中間値に基づき前記第１の粒子の密度を更新する
処理を含み、コンピュータにより実行される数値計算方法。

（付記６）
データ格納部と、
前記データ格納部に格納されているデータを用いて処理を行う処理部と、
を有し、
前記処理部は、
第１の基準密度の流体に係る第１の粒子の密度の代わりに、当該第１の粒子の密度の、前記第１の基準密度に対する比率をリーマン不変量において用いることで得られる第１の物理量を算出して、前記データ格納部に格納し、
第２の基準密度の流体に係る第２の粒子の密度の代わりに、当該第２の粒子の密度の、前記第２の基準密度に対する比率をリーマン不変量において用いることで得られる第２の物理量を算出して、前記データ格納部に格納し、
前記第１の物理量の勾配及び前記第２の物理量の勾配を算出して、前記データ格納部に格納し、
前記第１の粒子と前記第２の粒子との間における前記第１及び第２の物理量に関する時空間中間値を、前記第１の物理量、前記第２の物理量、前記第１の物理量の勾配及び前記第２の物理量の勾配を用いて算出して、前記データ格納部に格納し、
前記第１の粒子と前記第２の粒子との間における前記第１及び第２の物理量に関する中間値から、前記第１の粒子と前記第２の粒子との間における圧力の時空間中間値を算出して、前記データ格納部に格納し、
前記第１の粒子と前記第２の粒子との間における前記第１及び第２の物理量に関する中間値から、前記第１の粒子と前記第２の粒子との間における速度の時空間中間値を算出して、前記データ格納部に格納する
第１の処理を、前記第１の粒子との位置関係から特定される第２の粒子の各々について実行し、
算出された前記圧力の時空間中間値を用いて前記第１の粒子の速度を更新して、前記データ格納部に格納し、
算出された前記速度の時空間中間値に基づき前記第１の粒子の密度を更新して、前記データ格納部に格納する
処理を実行する情報処理装置。

１００ 情報処理装置
１１０ 初期条件データ格納部
１２０ 処理部
１３０ データ格納部
１４０ 処理結果格納部
１５０ 出力部

Claims (6)

1. 第１の基準密度の流体に係る第１の粒子の密度の代わりに、当該第１の粒子の密度の、前記第１の基準密度に対する比率をリーマン不変量において用いることで得られる第１の物理量を算出し、第２の基準密度の流体に係る第２の粒子の密度の代わりに、当該第２の粒子の密度の、前記第２の基準密度に対する比率をリーマン不変量において用いることで得られる第２の物理量を算出し、前記第１の物理量の勾配及び前記第２の物理量の勾配を算出し、前記第１の粒子と前記第２の粒子との間における前記第１及び第２の物理量に関する時空間中間値を、前記第１の物理量、前記第２の物理量、前記第１の物理量の勾配及び前記第２の物理量の勾配を用いて算出し、前記第１の粒子と前記第２の粒子との間における前記第１及び第２の物理量に関する時空間中間値から、前記第１の粒子と前記第２の粒子との間における圧力の時空間中間値を算出し、前記第１の粒子と前記第２の粒子との間における前記第１及び第２の物理量に関する時空間中間値から、前記第１の粒子と前記第２の粒子との間における速度の時空間中間値を算出する第１の処理を、前記第１の粒子との位置関係から特定される第２の粒子の各々について実行し、
   算出された前記圧力の時空間中間値を用いて前記第１の粒子の速度を更新し、
   算出された前記速度の時空間中間値に基づき前記第１の粒子の密度を更新する
   処理を、コンピュータに実行させるための数値計算プログラム。
2. 前記第１の物理量の勾配及び前記第２の物理量の勾配を算出する処理が、
   前記第１の基準密度の空間勾配を用いて、前記第１の物理量の勾配を算出し、
   前記第２の基準密度の空間勾配を用いて、前記第２の物理量の勾配を算出する
   処理を含む請求項１記載の数値計算プログラム。
3. 前記第１の粒子の速度を更新する処理が、
   前記圧力の時空間中間値から、前記第１の粒子の加速度を更新し、
   前記第１の粒子の加速度を用いて、前記第１の粒子の速度を更新する
   処理を含む請求項１又は２記載の数値計算プログラム。
4. 前記第１の粒子の密度を更新する処理が、
   前記第１の粒子の速度と前記速度の時空間中間値とを用いて、前記第１の粒子の密度時間変化を算出し、
   前記第１の粒子の密度時間変化を用いて、前記第１の粒子の密度を更新する
   処理を含む請求項１乃至３のいずれか１つ記載の数値計算プログラム。
5. 第１の基準密度の流体に係る第１の粒子の密度の代わりに、当該第１の粒子の密度の、前記第１の基準密度に対する比率をリーマン不変量において用いることで得られる第１の物理量を算出し、第２の基準密度の流体に係る第２の粒子の密度の代わりに、当該第２の粒子の密度の、前記第２の基準密度に対する比率をリーマン不変量において用いることで得られる第２の物理量を算出し、前記第１の物理量の勾配及び前記第２の物理量の勾配を算出し、前記第１の粒子と前記第２の粒子との間における前記第１及び第２の物理量に関する時空間中間値を、前記第１の物理量、前記第２の物理量、前記第１の物理量の勾配及び前記第２の物理量の勾配を用いて算出し、前記第１の粒子と前記第２の粒子との間における前記第１及び第２の物理量に関する時空間中間値から、前記第１の粒子と前記第２の粒子との間における圧力の時空間中間値を算出し、前記第１の粒子と前記第２の粒子との間における前記第１及び第２の物理量に関する時空間中間値から、前記第１の粒子と前記第２の粒子との間における速度の時空間中間値を算出する第１の処理を、前記第１の粒子との位置関係から特定される第２の粒子の各々について実行し、
   算出された前記圧力の時空間中間値を用いて前記第１の粒子の速度を更新し、
   算出された前記速度の時空間中間値に基づき前記第１の粒子の密度を更新する
   処理を含み、コンピュータにより実行される数値計算方法。
6. データ格納部と、
   前記データ格納部に格納されているデータを用いて処理を行う処理部と、
   を有し、
   前記処理部は、
   第１の基準密度の流体に係る第１の粒子の密度の代わりに、当該第１の粒子の密度の、前記第１の基準密度に対する比率をリーマン不変量において用いることで得られる第１の物理量を算出して、前記データ格納部に格納し、第２の基準密度の流体に係る第２の粒子の密度の代わりに、当該第２の粒子の密度の、前記第２の基準密度に対する比率をリーマン不変量において用いることで得られる第２の物理量を算出して、前記データ格納部に格納し、前記第１の物理量の勾配及び前記第２の物理量の勾配を算出して、前記データ格納部に格納し、前記第１の粒子と前記第２の粒子との間における前記第１及び第２の物理量に関する時空間中間値を、前記第１の物理量、前記第２の物理量、前記第１の物理量の勾配及び前記第２の物理量の勾配を用いて算出して、前記データ格納部に格納し、前記第１の粒子と前記第２の粒子との間における前記第１及び第２の物理量に関する時空間中間値から、前記第１の粒子と前記第２の粒子との間における圧力の時空間中間値を算出して、前記データ格納部に格納し、前記第１の粒子と前記第２の粒子との間における前記第１及び第２の物理量に関する時空間中間値から、前記第１の粒子と前記第２の粒子との間における速度の時空間中間値を算出して、前記データ格納部に格納する第１の処理を、前記第１の粒子との位置関係から特定される第２の粒子の各々について実行し、
   算出された前記圧力の時空間中間値を用いて前記第１の粒子の速度を更新して、前記データ格納部に格納し、
   算出された前記速度の時空間中間値に基づき前記第１の粒子の密度を更新して、前記データ格納部に格納する
   処理を実行する情報処理装置。

Images