(57)【要約】【課題】基準密度が異なる流体の混流現象を精度良く数値計算する。【解決手段】本方法は、第１基準密度の流体に係る第１粒子の密度の代わりに、当該第１粒子の密度の、第１基準密度に対する比率をリーマン不変量において用いることで得られる第１物理量を算出し、第２基準密度の流体に係る第２粒子の密度の代わりに、当該第２粒子の密度の、第２基準密度に対する比率をリーマン不変量において用いることで得られる第２物理量を算出し、第１及び第２物理量に基づき、第１粒子と第２粒子との間における第１及び第２物理量に関する時空間中間値を算出し、この時空間中間値から、第１粒子と第２粒子との間における圧力及び速度の時空間中間値を算出する処理を、第１粒子との位置関係から特定される第２粒子の各々について実行すれる。そして、圧力の時空間中間値を用いて第１粒子の速度を更新し、速度の時空間中間値に基づき第１粒子の密度を更新する。【選択図】図３
2014241002.tif 000002
【特許請求の範囲】
【請求項１】
第１の基準密度の流体に係る第１の粒子の密度の代わりに、当該第１の粒子の密度の、前記第１の基準密度に対する比率をリーマン不変量において用いることで得られる第１の物理量を算出し、
第２の基準密度の流体に係る第２の粒子の密度の代わりに、当該第２の粒子の密度の、前記第２の基準密度に対する比率をリーマン不変量において用いることで得られる第２の物理量を算出し、
前記第１の物理量の勾配及び前記第２の物理量の勾配を算出し、
前記第１の粒子と前記第２の粒子との間における前記第１及び第２の物理量に関する時空間中間値を、前記第１の物理量、前記第２の物理量、前記第１の物理量の勾配及び前記第２の物理量の勾配を用いて算出し、
前記第１の粒子と前記第２の粒子との間における前記第１及び第２の物理量に関する中間値から、前記第１の粒子と前記第２の粒子との間における圧力の時空間中間値を算出し、
前記第１の粒子と前記第２の粒子との間における前記第１及び第２の物理量に関する中間値から、前記第１の粒子と前記第２の粒子との間における速度の時空間中間値を算出する
第１の処理を、前記第１の粒子との位置関係から特定される第２の粒子の各々について実行し、
算出された前記圧力の時空間中間値を用いて前記第１の粒子の速度を更新し、
算出された前記速度の時空間中間値に基づき前記第１の粒子の密度を更新する
処理を、コンピュータに実行させるための数値計算プログラム。
【請求項２】
前記勾配を算出する処理が、
前記第１の基準密度の空間勾配を用いて、前記第１の物理量の勾配を算出し、
前記第２の基準密度の空間勾配を用いて、前記第２の物理量の勾配を算出する
処理を含む請求項１記載の数値計算プログラム。
【請求項３】
前記第１の粒子の速度を更新する処理が、
前記圧力の時空間中間値から、前記第１の粒子の加速度を更新し、
前記第１の粒子の加速度を用いて、前記第１の粒子の速度を更新する
処理を含む請求項１又は２記載の数値計算プログラム。
【請求項４】
前記第１の粒子の密度を算出する処理が、
前記第１の粒子の速度と前記速度の時空間中間値とを用いて、前記第１の粒子の密度時間変化を算出し、
前記第１の粒子の密度時間変化を用いて、前記第１の粒子の密度を更新する
処理を含む請求項１乃至３のいずれか１つ記載の数値計算プログラム。
【請求項５】
第１の基準密度の流体に係る第１の粒子の密度の代わりに、当該第１の粒子の密度の、前記第１の基準密度に対する比率をリーマン不変量において用いることで得られる第１の物理量を算出し、
第２の基準密度の流体に係る第２の粒子の密度の代わりに、当該第２の粒子の密度の、前記第２の基準密度に対する比率をリーマン不変量において用いることで得られる第２の物理量を算出し、
前記第１の物理量の勾配及び前記第２の物理量の勾配を算出し、
前記第１の粒子と前記第２の粒子との間における前記第１及び第２の物理量に関する時空間中間値を、前記第１の物理量、前記第２の物理量、前記第１の物理量の勾配及び前記第２の物理量の勾配を用いて算出し、
前記第１の粒子と前記第２の粒子との間における前記第１及び第２の物理量に関する中間値から、前記第１の粒子と前記第２の粒子との間における圧力の時空間中間値を算出し、
前記第１の粒子と前記第２の粒子との間における前記第１及び第２の物理量に関する中間値から、前記第１の粒子と前記第２の粒子との間における速度の時空間中間値を算出する
第１の処理を、前記第１の粒子との位置関係から特定される第２の粒子の各々について実行し、
算出された前記圧力の時空間中間値を用いて前記第１の粒子の速度を更新し、
算出された前記速度の時空間中間値に基づき前記第１の粒子の密度を更新する
処理を含み、コンピュータにより実行される数値計算方法。
【請求項６】
データ格納部と、
前記データ格納部に格納されているデータを用いて処理を行う処理部と、
を有し、
前記処理部は、
第１の基準密度の流体に係る第１の粒子の密度の代わりに、当該第１の粒子の密度の、前記第１の基準密度に対する比率をリーマン不変量において用いることで得られる第１の物理量を算出して、前記データ格納部に格納し、
第２の基準密度の流体に係る第２の粒子の密度の代わりに、当該第２の粒子の密度の、前記第２の基準密度に対する比率をリーマン不変量において用いることで得られる第２の物理量を算出して、前記データ格納部に格納し、
前記第１の物理量の勾配及び前記第２の物理量の勾配を算出して、前記データ格納部に格納し、
前記第１の粒子と前記第２の粒子との間における前記第１及び第２の物理量に関する時空間中間値を、前記第１の物理量、前記第２の物理量、前記第１の物理量の勾配及び前記第２の物理量の勾配を用いて算出して、前記データ格納部に格納し、
前記第１の粒子と前記第２の粒子との間における前記第１及び第２の物理量に関する中間値から、前記第１の粒子と前記第２の粒子との間における圧力の時空間中間値を算出して、前記データ格納部に格納し、
前記第１の粒子と前記第２の粒子との間における前記第１及び第２の物理量に関する中間値から、前記第１の粒子と前記第２の粒子との間における速度の時空間中間値を算出して、前記データ格納部に格納する
第１の処理を、前記第１の粒子との位置関係から特定される第２の粒子の各々について実行し、
算出された前記圧力の時空間中間値を用いて前記第１の粒子の速度を更新して、前記データ格納部に格納し、
算出された前記速度の時空間中間値に基づき前記第１の粒子の密度を更新して、前記データ格納部に格納する
処理を実行する情報処理装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、数値計算プログラム、数値計算方法及び情報処理装置に関する。
【背景技術】
【０００２】
流体や弾性体等の連続体の振る舞いを解く数値計算手法としては、格子をベースにして微分方程式の近似解を求める有限差分法、有限要素法、有限体積法等が知られている。近年では数値計算をＣＡＥ（Computer Aided Engineering）などの分野で活用するため、これらの数値計算手法も発展し、流体と構造物が相互作用する問題が解かれるようになってきた。しかしながら、格子を用いる手法では、自由表面などの界面の存在する問題や、流体−構造連成問題などの移動境界が発生する場合には取り扱いが複雑なため、プログラム作成が困難である場合が多い。
【０００３】
これに対して、格子を用いない粒子法（ＭＰＳ（Moving Particle Semi-implicit）法、ＳＰＨ（Smoothed Particle Hydrodynamics）法等）では、移動境界の取り扱いに特別な処置を必要としない。そのため近年、粒子法は広く用いられるようになってきている。
【０００４】
また、粒子法により高精度に計算する手法として、リーマンソルバを用いた計算手法が知られている。ある文献には、粒子法の計算で問題となる圧力の振動を抑制することで、精度よくシミュレーションを行う技術が開示されている。
【０００５】
例えば、１次元の非粘性流体の方程式におけるリーマン問題解を利用する計算手法（リーマンソルバと呼ぶ）を導入することにより、計算の高精度化を行うことを開示した文献が存在する。
【０００６】
ここで、リーマンソルバを粒子法に導入する方法の一例について簡単に説明する。僅かに圧縮する流体の方程式（０−１）乃至（０−３）を考える。
【数１】
2014241002.tif 000003
【０００７】
式（０−１）は質量保存則、式（０−２）が運動量保存則、式（０−３）は状態方程式である。ここでｖ（太字ｖ）は流体の速度ベクトル、ρは密度、ｐは圧力、ρ0は圧力が０となる密度（すなわち基準密度）である。また、ｃは音速であり、ｔは時刻を表す。
【０００８】
そして、式（０−１）乃至（０−３）を以下のように離散化する。
【数２】
2014241002.tif 000004
【０００９】
ここで、ｘ（太字ｘ）は３次元空間の位置ベクトルを表す。また、ｘi、ｖi、ρi及びｍiは、それぞれ、粒子ｉの位置ベクトル、速度ベクトル、密度、質量を表す。ｘijn+1/2は、粒子ｉと粒子ｊの相対位置ベクトルである。具体的には、以下のように表される。
ｘijn+1/2＝ｘin+1/2−ｘjn+1/2
【００１０】
また、上付き添え字はシミュレーションステップでの時刻を表わす。また、式（０−４）乃至（０−７）に含まれる記号は以下のように定義される。
【数３】
2014241002.tif 000005
【００１１】
さらに、Ｗはカーネル関数（重み関数とも呼ぶ）であり、以下のスプライン関数が良く用いられる。
【数４】
2014241002.tif 000006
【００１２】
ここでｈは、粒子間の影響半径で、初期状態の平均粒子間隔の２倍から３倍程度が良く用いられる。βはカーネル関数の全空間積分量が１になるように調整された値で、２次元の場合は０．７πｈ2 、３次元の場合はπｈ3と決められる。
【００１３】
式（０−８）及び（０−９）は、粒子ｉｊ間の中間値であって、１次元の流体の方程式を数値的に解き、ｄｔ／２だけ進んだ時刻における解を用いる。すなわち、図１に示すように、粒子ｉｊ間に軸ｌijをセットして、この軸上のベクトルに射影した速度ベクトルは、それぞれ以下のように表される。
【数５】
2014241002.tif 000007
【００１４】
具体的な決定方法は以下のようなものである。まず、式（０−１）乃至（０−３）について、空間次元が１次元の場合の方程式を整理する。
【数６】
2014241002.tif 000008
【００１５】
そして、式（０−１１）についてリーマン不変量ｑ+及びｑ-を用いて特性方程式の形式に書き直す。
【数７】
2014241002.tif 000009
【数８】
2014241002.tif 000010
【００１６】
そして、時刻ｎにおける粒子ｉ及びｊについてｑ+及びｑ-を以下のように表す。
【数９】
2014241002.tif 000011
【００１７】
さらに、空間勾配を以下のように求めておく。
【数１０】
2014241002.tif 000012
【００１８】
このようにして変数を整えた後に、図２に示すように、粒子ｉｊ間に物理量が不連続な面が存在する状況を初期値とした問題（リーマン問題）を考え、式（０−１２）及び（０−１３）の解析解を用いてｑ+及びｑ-の粒子ｉｊ間の中間における、時刻ｎ＋１／２の値を予測する。具体的には以下のように算出される。
【数１１】
2014241002.tif 000013
【００１９】
このような式（０−１６）及び（０−１７）を、式（０−５）及び（０−６）に代入することで粒子法のシミュレーションを高精度に実行できる。
【００２０】
一方、剛体と流体の連成計算や、混ざることのない別々の液体（たとえば水と油）の混流現象のシミュレーションは、自由表面や移動境界の扱いが簡単な粒子法が期待されている。混流現象の場合は安定状態となるときの密度が違う液体を取り扱うため、標準的なＳＰＨ法だけでは一般に計算は難しくなる。
【００２１】
密度の違う液体の混流を粒子法で扱う計算方法として、ある文献では、以下のような計算を行うことを示している。すなわち、密度の異なる流体方程式として以下の式（０−１８）乃至（０−２１）を用いる。
【数１２】
2014241002.tif 000014
【００２２】
ｇは重力による加速度であり、ρsは基準密度である。式（０−１）乃至（０−３）と異なる点は基準密度が位置によって異なることである。基準密度の空間変化で、異なる流体の混流を表している。
【００２３】
そして、式（０−１８）乃至（０−２１）をＳＰＨ法で以下のように各粒子について離散化する。
【数１３】
2014241002.tif 000015
【００２４】
式（０−２２）は質量保存則、式（０−２３）は運動量保存則、式（０−２４）は状態方程式を表す。ζ及びηはパラメータで定数である。ここでρ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.
【発明の概要】
【発明が解決しようとする課題】
【００２９】
従って、本発明の目的は、基準密度が異なる流体の混流現象を精度良く数値計算する技術を提供することである。
【課題を解決するための手段】
【００３０】
本発明に係る数値計算方法は、（Ａ）（ａ１）第１の基準密度の流体に係る第１の粒子の密度の代わりに、当該第１の粒子の密度の、第１の基準密度に対する比率をリーマン不変量において用いることで得られる第１の物理量を算出し、（ａ２）第２の基準密度の流体に係る第２の粒子の密度の代わりに、当該第２の粒子の密度の、第２の基準密度に対する比率をリーマン不変量において用いることで得られる第２の物理量を算出し、（ａ３）第１の物理量の勾配及び第２の物理量の勾配を算出し、（ａ４）第１の粒子と第２の粒子との間における第１及び第２の物理量に関する時空間中間値を、第１の物理量、第２の物理量、第１の物理量の勾配及び第２の物理量の勾配を用いて算出し、（ａ５）第１の粒子と第２の粒子との間における第１及び第２の物理量に関する中間値から、第１の粒子と第２の粒子との間における圧力の時空間中間値を算出し、（ａ６）第１の粒子と第２の粒子との間における第１及び第２の物理量に関する中間値から、第１の粒子と第２の粒子との間における速度の時空間中間値を算出する第１の処理を、第１の粒子との位置関係から特定される第２の粒子の各々について実行し、（Ｂ）算出された圧力の時空間中間値を用いて第１の粒子の速度を更新し、（Ｃ）算出された速度の時空間中間値に基づき第１の粒子の密度を更新する処理を含む。
【発明の効果】
【００３１】
一側面によれば、基準密度が異なる流体の混流現象を精度良く数値計算できるようになる。
【図面の簡単な説明】
【００３２】
【図１】図１は、粒子ｉｊの中間値を説明するための図である。
【図２】図２は、粒子ｉｊ間に物理量が不連続な面が存在する状況を表す図である。
【図３】図３は、本実施の形態において想定される混流の例を示す図である。
【図４】図４は、本実施の形態に係る情報処理装置の機能ブロック図である。
【図５】図５は、本実施の形態に係る処理の処理フローを示す図である。
【図６】図６は、本実施の形態に係る処理の処理フローを示す図である。
【図７】図７は、本実施の形態に係る処理の処理フローを示す図である。
【図８】図８は、本実施の形態に係る処理の処理フローを示す図である。
【図９】図９は、本実施の形態に係る処理の処理フローを示す図である。
【図１０】図１０は、本実施の形態に係る処理の処理フローを示す図である。
【図１１】図１１は、コンピュータの機能ブロック図である。
【発明を実施するための形態】
【００３３】
まず、本発明の実施の形態に係る演算の考え方について説明する。
【００３４】
本実施の形態では、以下の流体運動を表す方程式系を粒子法で解く。
【数１４】
2014241002.tif 000016
【００３５】
ここでｖ及びρは流体の速度場及び圧力場を表す。また、ρ及びρsは流体の密度及び基準密度（圧力が０の状態となる密度）を表す。また、ｃは音速を表す。
【００３６】
式（１）乃至（４）は、図３に示すような、水と油のように互いが混じり合うことなく、密度が異なる流体の混流を表す場合などに用いられる。
【００３７】
この式（１）乃至（４）を粒子法で空間離散化すると、以下のように表される。
【数１５】
2014241002.tif 000017
なお、ｊは例えば影響範囲内の他の粒子を表す。
【００３８】
さらに、時間について離散化することで以下の漸化式を導出する。
【数１６】
2014241002.tif 000018
【００３９】
ここで、ｄｔは時間刻み、その他の記号は以下のような意味である。
【数１７】
2014241002.tif 000019
【００４０】
式（１０）のように、密度の時間発展の右辺第二項におけるρni／ρnjの部分以外は、時刻ｎ＋１／２の値と時刻ｎとｎ＋１の平均値とを含むようになる。そうするとρni／ρnjの時間変化が小さい状況では、２次の精度を持つような計算方法となる。
【００４１】
ｖijlij,n+1/2及びｐijlij,n+1/2は粒子ｉｊ間の中間値であり、以下の手順で算出する。
【００４２】
まず、式（１）乃至（４）について１次元の方程式を考える。
【数１８】
2014241002.tif 000020
【００４３】
ここで、リーマン不変量におけるρに代わってρ／ρsを用いてリーマン不変量を変形することで、以下のような物理量を新たに導入する。
【数１９】
2014241002.tif 000021
【００４４】
そして、これらの物理量を用いて整理すると、以下のような式が得られる。
【数２０】
2014241002.tif 000022
【００４５】
式（１１−１）及び（１１−２）の右辺第一項だけを考えると、式（０−１２）及び（０−１３）と同様の形式であるため、右辺の第一項のみを考えて時間発展を以下のように行い、右辺第二項及び第三項を補う項を後で追加するようにして計算を行う。
Ａ）より具体的には、変数ｑ+及びｑ-を以下のようにして求める。
【数２１】
2014241002.tif 000023
これらは、時刻ｎにおける粒子ｉ及びｊについての上記物理量である。
Ｂ）また、変数ｑ+及びｑ-の勾配を以下のように計算する。
【数２２】
2014241002.tif 000024
【００４６】
式（１６）乃至（１９）における記号は以下のように表される。
【数２３】
2014241002.tif 000025
（２１）式は、変数ｑ+及びｑ-の勾配が、他の粒子の基準密度との関係に基づく基準密度の空間勾配を用いて算出されることを表している。すなわち、この分が、背景技術で述べた混流ではない場合との差である。
Ｃ）変数ｑ+及びｑ-について、粒子ｉｊ間の中間値を以下のように計算する。
【数２４】
2014241002.tif 000026
【００４７】
式（２３）及び（２４）の右辺第二項はそれぞれ式（０−１４）及び（０−１５）の右辺第二項と類似の形式であり、リーマンソルバの考え方を用いて得られた項であることが分かる。
【００４８】
なお、ρijn及びρs,ijnは粒子ｉｊ間での密度、及び基準密度の中間値で、以下のように表される。
【数２５】
2014241002.tif 000027
【００４９】
式（２３）及び（２４）の右辺第三項及び第四項は、それぞれ式（１１−１）及び（１１−２）の右辺第二項及び第三項を差分近似したものである。式（２３）式（２４）の右辺第三項のDρs /Dt|i およびDρs /Dt|j は、０として問題ない。また、基準密度の空間微分は、以下のように表され、標準的なＳＰＨ法の計算手法により得られる。
【数２６】
2014241002.tif 000028
【００５０】
Ｄ）上で算出した中間値ｑijn+1/2,+及びｑijn+1/2,-を用いて式（２７）及び（２８）より速度及び圧力の中間値を算出して、式（９）及び（１０）に代入する。
【数２７】
2014241002.tif 000029
【００５１】
以上の考え方に基づき、具体的な処理及び当該処理を行うための装置について説明する。
【００５２】
図４に、本実施の形態に係る情報処理装置１００の機能ブロック図を示す。情報処理装置１００は、初期条件データ格納部１１０と、処理部１２０と、データ格納部１３０と、処理結果格納部１４０と、出力部１５０とを有する。そして、情報処理装置１００は、出力装置２００に接続されている。
【００５３】
初期条件データ格納部１１０には、初期位置、初期速度、初期密度、基準密度、時間刻みｄｔ、粒子にかかる外力（例えば重力）のデータなどが格納されている。
【００５４】
処理部１２０は、初期条件データ格納部１１０に格納されているデータを用いて、各タイムステップについて、各粒子についての速度、密度及び位置等のデータを算出し、処理結果格納部１４０に格納する。なお、処理途中のデータについては、データ格納部１３０に格納される。
【００５５】
出力部１５０は、処理結果格納部１４０に格納されるデータを、出力装置２００（例えば印刷装置、表示装置、場合によってはネットワークで接続される他のコンピュータなどの装置）に出力する。
【００５６】
次に、図５乃至図１０を用いて、情報処理装置１００の処理内容について説明する。
【００５７】
まず、処理部１２０は、初期条件データ格納部１１０に格納されている初期条件データのデータを読み出す（図５：ステップＳ１）。また、処理部１２０は、ｎ＝１と設定する（ステップＳ３）。
【００５８】
そして、処理部１２０は、初期条件データにおいて設定がなされている粒子のうち未処理の粒子ｉを１つ特定する（ステップＳ５）。
【００５９】
その後、処理部１２０は、粒子ｉの位置をｄｔ／２だけ更新する（ステップＳ７）。具体的には、以下のような演算を実行する。ｘi1については、初期条件データに含まれる。
【数２８】
2014241002.tif 000030
【００６０】
また、処理部１２０は、粒子ｉの加速度、密度の時間変化項を初期化する（ステップＳ９）。具体的には、以下のように設定を行う（ステップＳ９）。
【数２９】
2014241002.tif 000031
【００６１】
さらに、処理部１２０は、粒子ｉにかかる外力を算出する（ステップＳ１１）。具体的には、以下のような演算を実行する。外力ｆniは、初期条件データに含まれる。
【数３０】
2014241002.tif 000032
【００６２】
また、処理部１２０は、粒子ｉの影響範囲に入る他の粒子ｊを抽出する（ステップＳ１３）。そして処理は端子Ａを介して図６の処理に移行する。
【００６３】
端子Ａを介して図６の処理の説明に移行して、処理部１２０は、抽出された粒子ｊの中から未処理の粒子ｊを１つ特定する（ステップＳ１５）。そして、処理部１２０は、粒子ｉｊ間の圧力に関する時空間中間値の算出処理を実行する（ステップＳ１７）。この処理については、図７及び図８を用いて説明する。
【００６４】
処理部１２０は、物理量ｑの中間値算出処理を実行する（図７：ステップＳ３１）。具体的には、図８の処理を実行する。
【００６５】
処理部１２０は、粒子ｉｊの密度、基準密度及び速度から、式（１３）及び（１４）に従って、物理量ｑ-i、ｑ+jを算出する（ステップＳ３５）。
【００６６】
また、処理部１２０は、粒子ｉｊの密度、基準密度及び速度の勾配から、式（１７）及び（１８）に従って、物理量ｑの勾配を算出する（ステップＳ３７）。
【００６７】
さらに、処理部１２０は、物理量ｑの時空間中間値ｑijn+1/2,+及びｑijn+1/2,-を、式（２３）及び（２４）に従って算出する（ステップＳ３９）。
【００６８】
図８の演算結果については、データ格納部１３０の容量が十分にあれば保持しておくことで、後の処理を省略できる。データ格納部１３０の容量が十分なければ、後の処理でも同じ演算を実行する。
【００６９】
図７の処理の説明に戻って、処理部１２０は、物理量ｑの時空間中間値を用いて、式（２８）に従って、圧力に関する時空間中間値ｐijn+1/2を算出する（ステップＳ３３）。
【００７０】
圧力に関する時空間中間値ｐijn+1/2を速度に関する時空間中間値とは別に先行して算出するのは、速度ｖin+1が、速度に関する時空間中間値の計算に用いられるためである。
【００７１】
図６の処理の説明に戻って、処理部１２０は、粒子ｉの加速度を、圧力に関する時空間中間値ｐijn+1/2を用いて更新する（ステップＳ１９）。このステップにおける演算は、以下のように表される。
【数３１】
2014241002.tif 000033
【００７２】
その後、処理部１２０は、ステップＳ１３で抽出された粒子のうち未処理の粒子ｊが存在しているか判断する（ステップＳ２１）。未処理の粒子ｊが存在する場合には、処理はステップＳ１５へ戻る。一方、未処理の粒子ｊが存在しない場合には、処理部１２０は、粒子ｉの加速度から粒子ｉの速度を更新する（ステップＳ２３）。このステップにおける演算は、以下のように表される。
【数３２】
2014241002.tif 000034
【００７３】
そして、処理部１２０は、ステップＳ１３で抽出された粒子ｊを未処理に設定する（ステップＳ２５）。なお、ステップＳ１３と同様の処理を実行するようにしても良い。処理は、端子Ｂを介して図９の処理に移行する。
【００７４】
図９の処理の説明に移行して、処理部１２０は、ステップＳ１３で抽出された粒子のうち未処理の粒子ｊを１つ特定する（ステップＳ４１）。そして、処理部１２０は、粒子ｉｊの速度に関する時空間中間値の算出処理を実行する（ステップＳ４３）。この処理については、図１０を用いて説明する。
【００７５】
処理部１２０は、物理量ｑの中間値算出処理を実行する（ステップＳ６３）。この処理は、図８の処理と同じである。上でも述べたように、物理量ｑの中間値については、データ格納部１３０に格納されている場合には、この処理を実行する代わりに、データ格納部１３０から読み出しても良い。
【００７６】
そして、処理部１２０は、物理量ｑの時空間中間値を用いて、式（２７）に従って、速度に関する時空間中間値ｖijlij,n+1/2を算出する（ステップＳ６５）。なお、式（２７）の後に以下の演算を行う。
【数３３】
2014241002.tif 000035
【００７７】
図９の処理の説明に戻って、処理部１２０は、粒子ｉの密度時間変化項を、以下の式に従って更新する（ステップＳ４５）。
【数３４】
2014241002.tif 000036
【００７８】
そして、処理部１２０は、ステップＳ１３で抽出された粒子のうち未処理の粒子ｊが存在するか判断する（ステップＳ４７）。未処理の粒子ｊが存在する場合には、処理はステップＳ４１に戻る。一方、未処理の粒子ｊが存在しない場合には、処理部１２０は、粒子ｉの密度時間変化項から粒子ｉの密度を、以下の式に従って更新する（ステップＳ４９）。
【数３５】
2014241002.tif 000037
【００７９】
さらに、処理部１２０は、粒子ｉの位置を、ステップＳ２３で算出された時刻ｎ＋１における速度ｖn+1iを用いてｄｔ／２だけ更新する（ステップＳ５１）。以下に示す式で算出する。
【数３６】
2014241002.tif 000038
【００８０】
そして、処理部１２０は、未処理の粒子ｉが存在するか判断する（ステップＳ５３）。未処理の粒子ｉが存在する場合には、端子Ｃを介して図５のステップＳ５に戻る。一方、未処理の粒子ｉが存在しない場合には、処理部１２０は、時刻ｎ＋１の結果を、処理結果格納部１４０に格納し、出力部１５０は、処理結果格納部１４０に格納された時刻ｎ＋１の処理結果を出力装置２００に対して所定の形式で出力する（ステップＳ５５）。
【００８１】
そして、処理部１２０は、時刻ｎ＋１は終了時刻であるか判断する（ステップＳ５７）。終了時刻でなければ、処理部１２０は、ｎ＝ｎ＋１と設定し（ステップＳ５９）、粒子ｉを全て未処理に設定する（ステップＳ６１）。そして、処理は端子Ｃを介して図５のステップＳ５に戻る。一方、時刻ｎ＋１が終了時刻であれば、処理を終了する。
【００８２】
以上のような処理を実行することで、粒子法で用いられているリーマンソルバによる計算の高精度化が実現する。
【００８３】
以上本発明の実施の形態を説明したが、本発明はこれに限定されるものではない。例えば、図４の機能ブロック図は一例であって、プログラムモジュール構成及びファイル構成とは一致しない場合もある。
【００８４】
処理フローについても、処理結果が変わらない限り、ステップの処理順番を入れ替えたり、複数ステップを並列に実行するようにしても良い。
【００８５】
さらに、情報処理装置１００は、１台のコンピュータではなく、複数のコンピュータで構成される場合もある。また、クライアントサーバ型システムで実装しても良いし、スタンドアロン型システムで実装しても良い。
【００８６】
なお、上で述べた情報処理装置１００は、コンピュータ装置であって、図１１に示すように、メモリ２５０１とＣＰＵ（Central Processing Unit）２５０３とハードディスク・ドライブ（ＨＤＤ：Hard Disk Drive）２５０５と表示装置２５０９に接続される表示制御部２５０７とリムーバブル・ディスク２５１１用のドライブ装置２５１３と入力装置２５１５とネットワークに接続するための通信制御部２５１７とがバス２５１９で接続されている。オペレーティング・システム（ＯＳ：Operating System）及び本実施例における処理を実施するためのアプリケーション・プログラムは、ＨＤＤ２５０５に格納されており、ＣＰＵ２５０３により実行される際にはＨＤＤ２５０５からメモリ２５０１に読み出される。ＣＰＵ２５０３は、アプリケーション・プログラムの処理内容に応じて表示制御部２５０７、通信制御部２５１７、ドライブ装置２５１３を制御して、所定の動作を行わせる。また、処理途中のデータについては、主としてメモリ２５０１に格納されるが、ＨＤＤ２５０５に格納されるようにしてもよい。本技術の実施例では、上で述べた処理を実施するためのアプリケーション・プログラムはコンピュータ読み取り可能なリムーバブル・ディスク２５１１に格納されて頒布され、ドライブ装置２５１３からＨＤＤ２５０５にインストールされる。インターネットなどのネットワーク及び通信制御部２５１７を経由して、ＨＤＤ２５０５にインストールされる場合もある。このようなコンピュータ装置は、上で述べたＣＰＵ２５０３、メモリ２５０１などのハードウエアとＯＳ及びアプリケーション・プログラムなどのプログラムとが有機的に協働することにより、上で述べたような各種機能を実現する。
【００８７】
以上述べた本実施の形態をまとめると、以下のようになる。
【００８８】
本実施の形態に係る数値計算方法は、（Ａ）（ａ１）第１の基準密度の流体に係る第１の粒子の密度の代わりに、当該第１の粒子の密度の、第１の基準密度に対する比率をリーマン不変量において用いることで得られる第１の物理量を算出し、（ａ２）第２の基準密度の流体に係る第２の粒子の密度の代わりに、当該第２の粒子の密度の、第２の基準密度に対する比率をリーマン不変量において用いることで得られる第２の物理量を算出し、（ａ３）第１の物理量の勾配及び第２の物理量の勾配を算出し、（ａ４）第１の粒子と第２の粒子との間における第１及び第２の物理量に関する時空間中間値を、第１の物理量、第２の物理量、第１の物理量の勾配及び第２の物理量の勾配を用いて算出し、（ａ５）第１の粒子と第２の粒子との間における第１及び第２の物理量に関する中間値から、第１の粒子と第２の粒子との間における圧力の時空間中間値を算出し、（ａ６）第１の粒子と第２の粒子との間における第１及び第２の物理量に関する中間値から、第１の粒子と第２の粒子との間における速度の時空間中間値を算出する第１の処理を、第１の粒子との位置関係から特定される第２の粒子の各々について実行し、（Ｂ）算出された圧力の時空間中間値を用いて第１の粒子の速度を更新し、（Ｃ）算出された速度の時空間中間値に基づき第１の粒子の密度を更新する処理を含む。
【００８９】
このようにリーマン不変量を変形して新たな物理量を定義することで、リーマンソルバによる高精度計算が実現される。
【００９０】
また、上で述べた勾配を算出する処理が、（ａ３１）第１の基準密度の空間勾配を用いて、第１の物理量の勾配を算出し、（ａ３２）第２の基準密度の空間勾配を用いて、第２の物理量の勾配を算出する処理を含むようにしても良い。リーマン不変量を、基準密度に対する粒子の密度の比率を用いて変形することで、リーマン不変量の勾配についても、第１の基準密度の空間勾配を用いるような形に変形される。
【００９１】
さらに、上で述べた第１の粒子の速度を更新する処理が、（ｂ１）圧力の時空間中間値から、第１の粒子の加速度を更新し、（ｂ２）第１の粒子の加速度を用いて、第１の粒子の速度を更新する処理を含むようにしても良い。加速度は、重力などの外力によって初期的に設定されるが、第２の粒子からの影響を重ね合わせて、さらに更新される。
【００９２】
また、上で述べた第１の粒子の密度を更新する処理が、（ｃ１）第１の粒子の速度と速度の時空間中間値とを用いて、第１の粒子の密度時間変化を算出し、（ｃ２）第１の粒子の密度時間変化を用いて、第１の粒子の密度を更新する処理を含むようにしても良い。
【００９３】
例えば、第１の処理の中で第１の粒子との位置関係から特定された第２の粒子について圧力に関する時空間中間値が得られて第１の粒子の速度が更新された後に、第１の粒子の密度時間変化を算出するようにしても良い。このようにすれば計算の精度が向上する。
【００９４】
なお、上で述べたような処理をコンピュータに実行させるためのプログラムを作成することができ、当該プログラムは、例えばフレキシブル・ディスク、ＣＤ−ＲＯＭなどの光ディスク、光磁気ディスク、半導体メモリ（例えばＲＯＭ）、ハードディスク等のコンピュータ読み取り可能な記憶媒体又は記憶装置に格納される。なお、処理途中のデータについては、ＲＡＭ等の記憶装置に一時保管される。
【００９５】
以上の実施例を含む実施形態に関し、さらに以下の付記を開示する。
【００９６】
（付記１）
第１の基準密度の流体に係る第１の粒子の密度の代わりに、当該第１の粒子の密度の、前記第１の基準密度に対する比率をリーマン不変量において用いることで得られる第１の物理量を算出し、
第２の基準密度の流体に係る第２の粒子の密度の代わりに、当該第２の粒子の密度の、前記第２の基準密度に対する比率をリーマン不変量において用いることで得られる第２の物理量を算出し、
前記第１の物理量の勾配及び前記第２の物理量の勾配を算出し、
前記第１の粒子と前記第２の粒子との間における前記第１及び第２の物理量に関する時空間中間値を、前記第１の物理量、前記第２の物理量、前記第１の物理量の勾配及び前記第２の物理量の勾配を用いて算出し、
前記第１の粒子と前記第２の粒子との間における前記第１及び第２の物理量に関する中間値から、前記第１の粒子と前記第２の粒子との間における圧力の時空間中間値を算出し、
前記第１の粒子と前記第２の粒子との間における前記第１及び第２の物理量に関する中間値から、前記第１の粒子と前記第２の粒子との間における速度の時空間中間値を算出する
第１の処理を、前記第１の粒子との位置関係から特定される第２の粒子の各々について実行し、
算出された前記圧力の時空間中間値を用いて前記第１の粒子の速度を更新し、
算出された前記速度の時空間中間値に基づき前記第１の粒子の密度を更新する
処理を、コンピュータに実行させるための数値計算プログラム。
【００９７】
（付記２）
前記勾配を算出する処理が、
前記第１の基準密度の空間勾配を用いて、前記第１の物理量の勾配を算出し、
前記第２の基準密度の空間勾配を用いて、前記第２の物理量の勾配を算出する
処理を含む付記１記載の数値計算プログラム。
【００９８】
（付記３）
前記第１の粒子の速度を更新する処理が、
前記圧力の時空間中間値から、前記第１の粒子の加速度を更新し、
前記第１の粒子の加速度を用いて、前記第１の粒子の速度を更新する
処理を含む付記１又は２記載の数値計算プログラム。
【００９９】
（付記４）
前記第１の粒子の密度を算出する処理が、
前記第１の粒子の速度と前記速度の時空間中間値とを用いて、前記第１の粒子の密度時間変化を算出し、
前記第１の粒子の密度時間変化を用いて、前記第１の粒子の密度を更新する
処理を含む付記１乃至３のいずれか１つ記載の数値計算プログラム。
【０１００】
（付記５）
第１の基準密度の流体に係る第１の粒子の密度の代わりに、当該第１の粒子の密度の、前記第１の基準密度に対する比率をリーマン不変量において用いることで得られる第１の物理量を算出し、
第２の基準密度の流体に係る第２の粒子の密度の代わりに、当該第２の粒子の密度の、前記第２の基準密度に対する比率をリーマン不変量において用いることで得られる第２の物理量を算出し、
前記第１の物理量の勾配及び前記第２の物理量の勾配を算出し、
前記第１の粒子と前記第２の粒子との間における前記第１及び第２の物理量に関する時空間中間値を、前記第１の物理量、前記第２の物理量、前記第１の物理量の勾配及び前記第２の物理量の勾配を用いて算出し、
前記第１の粒子と前記第２の粒子との間における前記第１及び第２の物理量に関する中間値から、前記第１の粒子と前記第２の粒子との間における圧力の時空間中間値を算出し、
前記第１の粒子と前記第２の粒子との間における前記第１及び第２の物理量に関する中間値から、前記第１の粒子と前記第２の粒子との間における速度の時空間中間値を算出する
第１の処理を、前記第１の粒子との位置関係から特定される第２の粒子の各々について実行し、
算出された前記圧力の時空間中間値を用いて前記第１の粒子の速度を更新し、
算出された前記速度の時空間中間値に基づき前記第１の粒子の密度を更新する
処理を含み、コンピュータにより実行される数値計算方法。
【０１０１】
（付記６）
データ格納部と、
前記データ格納部に格納されているデータを用いて処理を行う処理部と、
を有し、
前記処理部は、
第１の基準密度の流体に係る第１の粒子の密度の代わりに、当該第１の粒子の密度の、前記第１の基準密度に対する比率をリーマン不変量において用いることで得られる第１の物理量を算出して、前記データ格納部に格納し、
第２の基準密度の流体に係る第２の粒子の密度の代わりに、当該第２の粒子の密度の、前記第２の基準密度に対する比率をリーマン不変量において用いることで得られる第２の物理量を算出して、前記データ格納部に格納し、
前記第１の物理量の勾配及び前記第２の物理量の勾配を算出して、前記データ格納部に格納し、
前記第１の粒子と前記第２の粒子との間における前記第１及び第２の物理量に関する時空間中間値を、前記第１の物理量、前記第２の物理量、前記第１の物理量の勾配及び前記第２の物理量の勾配を用いて算出して、前記データ格納部に格納し、
前記第１の粒子と前記第２の粒子との間における前記第１及び第２の物理量に関する中間値から、前記第１の粒子と前記第２の粒子との間における圧力の時空間中間値を算出して、前記データ格納部に格納し、
前記第１の粒子と前記第２の粒子との間における前記第１及び第２の物理量に関する中間値から、前記第１の粒子と前記第２の粒子との間における速度の時空間中間値を算出して、前記データ格納部に格納する
第１の処理を、前記第１の粒子との位置関係から特定される第２の粒子の各々について実行し、
算出された前記圧力の時空間中間値を用いて前記第１の粒子の速度を更新して、前記データ格納部に格納し、
算出された前記速度の時空間中間値に基づき前記第１の粒子の密度を更新して、前記データ格納部に格納する
処理を実行する情報処理装置。
【符号の説明】
【０１０２】
１００  情報処理装置
１１０  初期条件データ格納部
１２０  処理部
１３０  データ格納部
１４０  処理結果格納部
１５０  出力部
【図１】
2014241002.tif 000039
【図２】
2014241002.tif 000040
【図３】
2014241002.tif 000041
【図４】
2014241002.tif 000042
【図５】
2014241002.tif 000043
【図６】
2014241002.tif 000044
【図７】
2014241002.tif 000045
【図８】
2014241002.tif 000046
【図９】
2014241002.tif 000047
【図１０】
2014241002.tif 000048
【図１１】
2014241002.tif 000049