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

Description

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

数値計算を用いて水や空気等の流体の流れを調べる流体解析や、圧縮されたゴム等の弾性体の振る舞いを調べる弾性体解析において、図１に模式的に示すように、連続体を粒子の分布で表現する手法（以下では総称して粒子法と呼ぶ）が提案されている。そして、粒子法の代表的なものとして、Smoothed Particles Hydrodynamics（ＳＰＨ）法や、Moving Particles Semi-implicit（ＭＰＳ）法が知られている。粒子法における標準的な手法としては、図２に示すように、ある粒子ａに対して予め領域（以下では、「影響領域」と呼ぶ。図２の例では半径２ｈ（ｈは粒子の影響範囲の大きさを表すパラメータ（以下、影響半径と呼ぶ））で規定される円。）を設定しておき、その影響領域の内に存在する相手粒子ｂのみからの相互作用を計算することで物体の運動を解析する。

例えば、鋳造、鍛造などの金属加工では、液体金属の中に冷えて固まった金属が混じる固化過程において体積変化を伴うなど、複雑な過程となる。より具体的には、図３に示すように、容器内部に液体金属を充填すると、容器内部の液体金属の熱は、熱伝導で外側方向に伝わると共に、容器表面及び上面の液体金属表面では放射冷却が発生する。このような現象が進むと、液体金属が固化することになる。図３右側の例では、固化による体積縮小で一部に空間が生じてしまう。

自由表面の扱いの簡単さや複数の計算ノードで並列計算を行う場合の並列性能面、固体との連成計算が比較的容易であるといった利点により鋳造や鍛造シミュレーションは粒子法の活躍が期待される分野であると考えられる。

上で述べたＭＰＳ法については、樹脂流動解析に用いた例が存在している。

また、上で述べたＳＰＨ法を用いた例も存在しており、この例では、複数の粒子の物理量をカーネル関数と呼ぶ重み関数を用いて平滑化して、基礎方程式を離散化している。例えば、鋳造における溶融金属の流れを扱うために、連続を表す式、運動方程式、及びエネルギー方程式を、以下のように離散化している。

ここで、下付け添え字は、各粒子のインデックスを表す。ρ_aは粒子ａの密度を表し、ｍ_aは粒子ａの質量を表し、ｖ_aは粒子ａの速度ベクトルを表し、ｒ_aは粒子ａの位置ベクトルを表し、Ｐ_aは粒子ａの圧力を表し、μ_aは粒子ａの粘性係数を表し、ｕ_aは粒子ａの単位質量あたりの内部エネルギーを表し、Ｔ_aは粒子ａの温度を表し、κ_aは粒子ａの熱伝導係数を表す。

さらに、以下の関係が成り立つ。（図４参照のこと）
ｖ_ab＝ｖ_a−ｖ_b
ｒ_ab（太字）＝ｒ_a−ｒ_b
ｒ_ab（イタリック）＝｜ｒ_ab（太字）｜
Ｔ_ab＝Ｔ_a−Ｔ_b

また、ηは分母の発散を抑える数値パラメータを表し、ｇは重力加速度を表し、ξは粘性の効果を示す数値パラメータを表し、ξ＝4.96333を用いると良いとされている。Ｗはカーネル関数で、粒子の分布から連続場を構成するのに用いる。例えば、Ｗは、３次のスプライン関数が用いられる。

また、圧力については、以下のような状態方程式を用いて求める。

但し、Ｐ₀はその流体の基準圧力を表し、ρ₀は基準密度を表し、γは断熱指数を表す。γ＝１やγ＝７が良く用いられる。

エネルギー方程式である（３）式の右辺は、熱伝導に由来する項である。断熱膨張、圧縮及び粘性による散逸で運動エネルギーが内部エネルギーに転換される効果を考慮する場合には、以下の式を（３）式に加える。

鋳造などの過程では物体の温度変化を考える上で、図３のような例では、熱伝導と、断熱膨張、圧縮、粘性散逸の項に加えて物体表面からの熱放射による温度変化を考慮することになる。具体的には、黒体の表面から単位面積、単位時間あたりに放出される電磁波のエネルギーが、その黒体の熱力学温度Ｔの４乗に比例するという物理法則であるシュテファン・ボルツマンの法則を考慮に入れることになる。

しかしながら、放射冷却について、粒子の様々な状態を考慮した形で内部エネルギーを算出するような技術は存在していない。具体的には、例えば図５のような状況では、とがった形状の先端（丸印の部分）では、図６のように物体表面が平面となっている場合における当該平面上の粒子（表面粒子と呼ぶ。なお、表面より内側の粒子を非表面粒子又は内部粒子と呼ぶ）よりも放射量が多いはずであり、平面上に存在する粒子よりも単位時間あたりに失われるエネルギーは大きいはずである。しかし、このような違いを数値計算に反映させることが考慮されていない。特に、液滴が飛び散るなどして、孤立粒子になった場合（例えば図７）は全空間に向けて熱放射を行うため、放射冷却が効果的に起きるが、従来の数値計算においてはそのような効果が考慮されていない。

Paul W. Cleary,"Extension of SPH to predict feeding, freezing and defect creation in low pressure die casting", Applied Mathematical Modelling, vol.34, (2010), pp.3189-3201 Paul W. Cleary,"Modelling confined multi-material heat and mass flows using SPH", Applied Mathematical Modelling, vol.22, (1998), pp.981-993

特開２００２−１３７２７２号公報 特開２０１２−１５０６７３号公報 国際公開第２０１２／１１１０８２号パンフレット

本発明の目的は、一側面によれば、粒子の状態を考慮して放射冷却についての数値計算を可能にするための技術を提供することである。

本発明の一態様に係る数値計算方法は、（Ａ）連続体を粒子の集合として表現する場合における各粒子について、当該粒子が連続体表面にさらされている度合いに対応する係数に応じ且つ放射冷却に基づく内部エネルギーの時間微分を算出し、（Ｂ）各粒子について算出された内部エネルギーの時間微分に基づき、単位時間後の内部エネルギーを算出する処理を含む。

一側面によれば、粒子の状態を考慮して放射冷却についての数値計算を行うことが可能となる。

図１は、粒子の集合によって連続体を表す場合を示す図である。 図２は、影響範囲を説明するための図である。 図３は、固化過程の一例を示す図である。 図４は、パラメータの説明のための図である。 図５は、粒子の状態を表す図である。 図６は、粒子の状態を表す図である。 図７は、粒子の状態を表す図である。 図８は、粒子配置と、規格化された数密度と、多項式の変数との関係を表す図である。 図９は、Ｑ_aとｇ_a（Ｑ_a）の関係を表す図である。 図１０は、本実施の形態に係る情報処理装置の機能ブロック図である。 図１１は、本実施の形態に係る処理の処理フローを示す図である。 図１２は、本実施の形態に係る処理の処理フローを示す図である。 図１３は、コンピュータの機能ブロック図である。

最初に、本発明の実施の形態に係る考え方を述べる。

具体的には、全粒子について、以下の式に従って内部エネルギーの時間微分を算出する。

ここでｇ_aは、粒子ａが物体表面（連続体表面）にさらされている度合いを定量化するための係数であって、
（Ａ）粒子が連続体内部に存在する場合には放射冷却が起こらず、
（Ｂ）粒子が平面上に一様分布している場合には放射冷却が、当該粒子と等価な体積を有する、各辺が等長である平行六面体が有する六面のうち一面からの放射冷却が起こり、且つ
（Ｃ）粒子が孤立している場合には当該粒子と等価な体積の球と同じ放射冷却が起こる
ように補正するための連続関数である。なお、シュテファン・ボルツマンの法則は、（８）式において、係数ｇ_a以外の部分で表される。また、σは、物質に依存しないシュテファン・ボルツマン係数であり、αは物体の、黒体からのずれを表す係数（以下、放射率と呼ぶ）を表し、ｄは空間次元を表す。ｄは、１、２又は３である。さらに、（ｍ_a／ρ_a）^1-(1/d)は、粒子ａの表面積を表している。ここで、「当該粒子と等価な体積を有する、各辺が等長である平行六面体が有する六面のうち一面」は、当該粒子と等価な体積の立方体が有する六面のうち一面と等しい面積を有することになる。また、係数ｇ_aの満たすべき（Ｂ）の条件の代わりとして、表面形状に応じて粒子を、当該粒子と等価な体積の多面体や半球と仮定した場合において、当該粒子と等価な体積の立方体が有する六面のうち一面と等しい面積を有する立体であっても良い。

本実施の形態では、このため、係数ｇ_aを、以下に示すような、規格化された数密度の関数として表す。

この規格化された数密度は、着目している粒子ａの影響範囲内に他の粒子が十分にあると１に近い値となり、影響範囲内の粒子が少ない場合には小さな値となる。

ここで、カーネル関数Ｗ（ｒ，ｈ）は、例えば以下のような５次のスプライン関数である。

ここでα_dは規格化係数である。また、ｑ＝ｒ／ｈである。

規格化された数密度ｓ_aは、物体内部では「１」、孤立した粒子では以下の値となる。

さらに、平面とみなせる平面上では「（１＋ｓ_min,a）／２」という値になるように定義される。

このような条件を満たすような係数ｇ_aは一意には決まらないが、例えば、以下のような多項式を採用すれば、上で述べた（Ａ）乃至（Ｃ）の条件は満たされる。

ここで、記号は以下のような意味である。

以下、この式で（Ａ）乃至（Ｃ）の条件を満たしていることを説明する。

Ｑ_a＝（１−ｓ_a）／（１−ｓ_min,a）と定義したとき、粒子配置と、ｓ_a及びＱ_aの値の関係は、図８に示すようになる。すなわち、内部粒子（＝非表面粒子）の場合には、ｓ_a＝１で、Ｑ_a＝０である。また、表面粒子の場合には、ｓ_a＝（１＋ｓ_min,a）／２で、Ｑ_a＝１／２である。さらに、孤立粒子の場合には、ｓ_a＝ｓ_min,aで、Ｑ_a＝１である。

ここで係数ｇ_aの満たすべき条件を、多項式ｇ_a（Ｑ_a）の満たすべき条件として述べると、以下のようになる。
（Ａ）粒子が物体内部に存在するときには放射冷却は起こらない。
これはｇ_a（０）＝０と等価である。なぜならば、内部粒子の時Ｑ_a＝０であり、この時ｇ_a（０）＝０であれば、放射冷却によって失われる内部エネルギーが０となり、放射冷却が起きていないことになるからである。
（Ｂ）平面上に一様に分布しているとき、当該粒子と等価な体積の立方体のうち一面からの放射冷却が起こる（従来手法）。
これはｇ_a（１／２）＝１と等価である。なぜならば、表面粒子の場合には、Ｑ_a＝１／２であり、この時ｇ_a（１／２）＝１であれば、放射冷却によって失われる内部エネルギーは、シュテファン・ボルツマンの法則に従って計算した量（（８）式においてｇ_aを除外した式の値）に一致するためである。
（Ｃ）孤立した粒子は等価な体積の球と同様に放射冷却が生ずる。
この条件は、以下のような式で表される。

ここでＳ_1,aは、上で述べたとおりに定義され、粒子ａを球とみなしたときの表面積を表す。（ｍ_a／ρ_a）^1-(1/d)は平面上に一様に分布している配置での粒子ａ１つあたりの表面積を表す。孤立した粒子については、これに係数ｇ_aを掛けた結果が、粒子ａを球とみなした時の表面積と一致させたならば、等価な体積の球と同様に放射冷却をするという条件を満たすことになる。

ここで、数値計算上、内部粒子についてはＱ_aの値が粒子配置の多少のずれによって０付近で多少揺らぐと考えられるが、それによって放射冷却の補正係数が敏感に変動することを避けることが好ましい。従って、下記の式の条件を付加する。

よって、以下に示す４式を満たし且つ最も次数の低い多項式を求めると（９）式となる。

但し、最も次数の低い多項式ではなく、図９に示すような滑らかな関数となれば、より次数の高い多項式などを採用するようにしても良い。

なお、３次元においてはＡ_a＝３^2/3（４π）^1/3となるので、この時のＱ_aとｇ_a（Ｑ_a）との関係は、図９に示すようになる。

図９では、Ｑ_a＝０付近ではｇ_a（Ｑ_a）は滑らかに増加し、Ｑ_a＝１／２ではｇ_a（Ｑ_a）＝１となっている。Ｑ_aがそれより大きくなると、ｇ_a（Ｑ_a）は急激に増加する。

このような係数ｇ_aを採用することによって、粒子の状態を放射冷却による内部エネルギーの減少に適切に反映させることができる。また、このような係数ｇ_aは、規格化された数密度を変数とする多項式で表される連続関数であって、処理フローにおける分岐といった形ではなく無段階で粒子の状態を反映させることができるようになっている。

次に、このような放射冷却に係るシミュレーションを行うための情報処理装置１００の機能ブロック図を図１０に示す。

本実施の形態に係る情報処理装置１００は、入力部１１０と、第１データ格納部１２０と、物理量算出部１３０と、第２データ格納部１４０と、出力部１５０とを有する。

入力部１１０は、例えばネットワークで接続されている他のコンピュータからデータを取得して又はユーザからのデータ入力を受け付けて、処理すべきデータとして第１データ格納部１２０に格納する。

入力データは、数値計算を行う対象である連続体粒子のデータと、連続体粒子の運動に関する境界条件などを設定する固定境界要素のデータとを含む。連続体粒子は、流体をモデル化したものであって、これを表現するためのデータは、例えば初期の中心位置座標、初期の速度、影響半径、密度、質量、粘性等である。固定境界要素は、鋳型などの表面などを微小部分に分割した面要素をモデル化したものであって、このデータとしては、例えば微小な円盤の集合として境界全体をモデル化して、各境界要素に中心座標、面要素の法線ベクトル、面要素の面積を含む。ポリゴンの集合として境界全体を表現して、各境界要素に複数の頂点の位置座標を設定するようにしても良い。

物理量算出部１３０は、近傍リスト生成部１３１と、放射冷却算出部１３２と、積分処理部１３３とを有する。物理量算出部１３０は、単位時間毎に各粒子の物理量を算出する処理を行う。そのために、近傍リスト生成部１３１は、単位時間毎に、各粒子の影響範囲に入っている他の粒子のリストを生成し、第２データ格納部１４０に格納する。また、放射冷却算出部１３２は、単位時間毎に、各粒子について放射冷却による内部エネルギーの時間微分を算出し、第２データ格納部１４０に格納する。積分処理部１３３は、各粒子について算出された加速度、速度、密度の時間微分、内部エネルギーの時間微分等を時間積分して、１単位時間後の物理量（例えば速度、位置、密度、内部エネルギー等）を算出し、第２データ格納部１４０に格納する。

出力部１５０は、第２データ格納部１４０に格納されている単位時間毎の物理量を用いて出力データを生成し、他のコンピュータに出力したり、印刷装置や表示装置などの出力装置に出力する。

次に、図１１及び図１２を用いて情報処理装置１００の処理内容について説明する。まず、入力部１１０は、他のコンピュータから処理に用いられるデータを取得するか、ユーザからの入力を受け付けて、第１データ格納部１２０に格納する（ステップＳ１）。上で述べたような入力データを取得して、第１データ格納部１２０に格納する。

次に、物理量算出部１３０は、時刻ｔを０に初期化し（ステップＳ３）、近傍リスト生成部１３１は、時刻ｔにおける粒子分布に基づいて、各粒子について近傍粒子リストを生成し、第２データ格納部１４０に格納する（ステップＳ５）。

ｔ＝０の場合には、第１データ格納部１２０に格納されているデータに含まれる初期位置から、例えば着目している粒子から距離が影響範囲２ｈ以内の他の粒子の識別子を、リストに追加する。影響範囲は、相互に影響を与え合う範囲であり、例えば、粒子が移動する場合に、他の粒子へ力を加える等の処理を行う範囲である。なお、ｔ＝１以降については、第２データ格納部１４０に格納されているデータを用いて、近傍粒子リストを生成する。

そして、物理量算出部１３０は、未処理の粒子ａを１つ特定する（ステップＳ７）。その後、物理量算出部１３０は、連続体（物体）の物理モデル（例えば（１）式乃至（３）式）に応じて、特定された粒子ａについての解析対象量（放射冷却以外の要素に由来する内部エネルギーの時間微分を含む）を、当該粒子及びその近傍粒子リストに含まれる粒子との相互作用の重ね合わせとして算出し、第２データ格納部１４０に格納する（ステップＳ９）。このステップの処理については、例えば非特許文献２などで説明されており、ここでは説明を省略する。

解析対象量としては、例えば加速度、速度、密度の時間微分、放射冷却以外の要素に由来する内部エネルギーの時間微分が含まれる。放射冷却以外の要素に由来する内部エネルギーの時間微分は、以下のように表される。

ｕ_aは、粒子ａの内部エネルギーを表す。

さらに、放射冷却算出部１３２は、特定された粒子ａについて、物体表面にさらされている度合いに応じた補正係数ｇ_aを用いて、放射冷却によって単位時間あたりに失われる内部エネルギー量（内部エネルギーの時間微分）を算出し、第２データ格納部１４０に格納する（ステップＳ１１）。具体的には、（８）式及び（９）式に従って、内部エネルギーの時間微分を算出する。ここでは、放射冷却に由来する内部エネルギーの時間微分は、以下のように表される。

そして、物理量算出部１３０は、算出された内部エネルギー量を、内部エネルギーの時間微分に加算し、第２データ格納部１４０に格納する（ステップＳ１３）。具体的には、以下のような演算を行う。

そして、物理量算出部１３０は、未処理の粒子が存在するか判断する（ステップＳ１５）。未処理の粒子が存在する場合には、処理はステップＳ７に戻る。一方、未処理の粒子が存在しない場合には、処理は端子Ａを介して図１２の処理に移行する。

図１２の処理の説明に移行して、物理量算出部１３０の積分処理部１３３は、各粒子について、解析対象量の時間積分を実行し、処理結果を第２データ格納部１４０に格納する（ステップＳ１７）。着目している内部エネルギーについては、以下のような演算を実行する。

ここでｕ^t _aは時刻ｔにおける内部エネルギーを表している。

その他の解析対象量についても、同様の演算を行って、各粒子について、時刻ｔ＋１についての物理量を算出する。

そして、出力部１５０は、第２データ格納部１４０に格納されている時刻ｔ＋１についての物理量を、出力装置（他のコンピュータ、印刷装置又は表示装置等）に出力する（ステップＳ１９）。

物理量算出部１３０は、例えば時刻ｔが処理終了時刻になったか否かを判断することで処理終了であるか否かを判断する（ステップＳ２１）。処理終了であれば、処理を終了し、処理終了でなければ、物理量算出部１３０は、時刻ｔを１インクリメントして（ステップＳ２３）、処理は端子Ｂを介してステップＳ５に移行する。

以上のような処理を行えば、係数ｇ_aは、粒子毎にその状態を反映した値、すなわち、表面にさらされている度合いに応じた値となる。そして、各粒子が、単位時間毎に異なる状態であっても、その状態に応じた放射冷却による内部エネルギーの時間微分が算出できるようになる。

この際、処理フローの分岐という形ではなく、連続関数として補正係数ｇ_aが定義されているので、放射冷却の効果を無段階に反映させることができ、正確に物体の温度変化を計算することができるようになる。

以上本発明の実施の形態を説明したが、本発明はこれに限定されるものではない。例えば、図１０に示した機能ブロック図は一例であって、プログラムモジュール構成とは一致しない場合がある。また、処理フローについても、各粒子の演算については複数のプロセッサで並列演算を行うように処理フローを変形したりするようにしても良い。

また、情報処理装置１００は、１台のコンピュータではなく、複数台のコンピュータでその機能が分担される場合もある。

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

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

本実施の形態に係る数値計算方法は、（Ａ）連続体を粒子の集合として表現する場合における各粒子について、当該粒子が連続体表面にさらされている度合いに対応する係数に応じ且つ放射冷却に基づく内部エネルギーの時間微分を算出し、（Ｂ）各粒子について算出された内部エネルギーの時間微分に基づき、単位時間後の内部エネルギーを算出する処理を含む。

このようにすれば、個々の粒子の状態を反映させることができ、状態変化が生じてもその状態変化に応じて、放射冷却に基づき内部エネルギーの時間微分を算出できるようになる。

上で述べた係数が、（ａ）粒子が連続体内部に存在する場合には放射冷却が起こらず、（ｂ）粒子が平面上に一様分布している場合には当該粒子と等価な体積の立方体のうち一面からの放射冷却が起こり、且つ（ｃ）粒子が孤立している場合には当該粒子と等価な体積の球と同じ放射冷却が起こるように補正するための連続関数である場合もある。

このようにすれば放射冷却による内部エネルギーの時間微分を正確に算出できるようになる。

さらに、上で述べた連続関数が、規格化された数密度ｓ_aの関数となっている

（ｂは、粒子ａの所定範囲内の各粒子を表し、ｍは粒子の質量を表し、ρは密度を表し、ｒ_abは粒子ａと粒子ｂの距離を表し、ｈは粒子の影響範囲の大きさを表すパラメータであり、Ｗはカーネル関数を表す。）場合もある。規格化された数密度は、影響範囲内の粒子ｂに応じて変化するため、粒子ａの状態に応じた放射冷却による内部エネルギーの時間微分が算出されるようになる。

さらに、上で述べた連続関数が、

（ｓ_min,a＝ｍ_a／ρ_aＷ（０，ｈ））
の多項式で表される場合もある。このような変数変換を行うことで、多項式についての条件設定が容易になる。

さらに、上で述べた連続関数が、

である場合もある。上で述べた（ａ）乃至（ｃ）の３条件に対応する条件とＱ_a＝０付近で滑らかであるという条件とを満たす次数の最も小さいものを採用することで、滑らかに変化する連続関数を用いることができるようになる。

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

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

（付記１）
連続体を粒子の集合として表現する場合における各粒子について、当該粒子が前記連続体表面にさらされている度合いに対応する係数に応じ且つ放射冷却に基づく内部エネルギーの時間微分を算出し、
各粒子について算出された前記内部エネルギーの時間微分に基づき、単位時間後の前記内部エネルギーを算出する
処理を、コンピュータに実行させるための数値計算プログラム。

（付記２）
前記係数が、
粒子が前記連続体内部に存在する場合には放射冷却が起こらず、
粒子が平面上に一様分布している場合には当該粒子と等価な体積を有する、各辺が等長である平行六面体が有する六面のうち一面からの放射冷却が起こり、且つ
粒子が孤立している場合には当該粒子と等価な体積の球と同じ放射冷却が起こる
ように補正するための連続関数である
付記１記載の数値計算プログラム。

（付記３）
前記連続関数が、規格化された数密度ｓ_aの関数となっている

（ｂは、粒子ａの所定範囲内の各粒子を表し、ｍは粒子の質量を表し、ρは密度を表し、ｒ_abは粒子ａと粒子ｂの距離を表し、ｈは粒子の影響範囲の大きさを表すパラメータであり、Ｗはカーネル関数を表す。）
付記２記載の数値計算プログラム。

（付記４）
前記連続関数が、

（ｓ_min,a＝ｍ_a／ρ_aＷ（０，ｈ））
の多項式で表される付記３記載の数値計算プログラム。

（付記５）
前記連続関数が、

である付記３記載の数値計算プログラム。

（付記６）
連続体を粒子の集合として表現する場合における各粒子について、当該粒子が前記連続体表面にさらされている度合いに対応する係数に応じ且つ放射冷却に基づく内部エネルギーの時間微分を算出し、
各粒子について算出された前記内部エネルギーの時間微分に基づき、単位時間後の前記内部エネルギーを算出する
処理を含み、コンピュータにより実行される数値計算方法。

（付記７）
連続体を粒子の集合として表現する場合における各粒子について、当該粒子が前記連続体表面にさらされている度合いに対応する係数に応じ且つ放射冷却に基づく内部エネルギーの時間微分を算出する手段と、
各粒子について算出された前記内部エネルギーの時間微分に基づき、単位時間後の前記内部エネルギーを算出する手段と、
を有する情報処理装置。

１１０ 入力部
１２０ 第１データ格納部
１３０ 物理量算出部
１４０ 第２データ格納部
１５０ 出力部
１３１ 近傍リスト生成部
１３２ 放射冷却算出部
１３３ 積分処理部

Claims (6)

1. 連続体を粒子の集合として表現する場合における各粒子について、当該粒子が前記連続体表面にさらされている度合いに対応する係数に応じ且つ放射冷却に基づく内部エネルギーの時間微分を算出し、
   各粒子について算出された前記内部エネルギーの時間微分に基づき、単位時間後の前記内部エネルギーを算出する
   処理を、コンピュータに実行させるための数値計算プログラムであって、
   前記係数が、
   粒子が前記連続体内部に存在する場合には放射冷却が起こらず、
   粒子が平面上に一様分布している場合には当該粒子と等価な体積を有する、各辺が等長である平行六面体が有する六面のうち一面からの放射冷却が起こり、且つ
   粒子が孤立している場合には当該粒子と等価な体積の球と同じ放射冷却が起こる
   ように補正するための連続関数である
   数値計算プログラム。
2. 前記連続関数が、規格化された数密度ｓ_aの関数となっている
   

   

   （ｂは、粒子ａの所定範囲内の各粒子を表し、ｍは粒子の質量を表し、ρは密度を表し、ｒ_abは粒子ａと粒子ｂの距離を表し、ｈは粒子の影響範囲の大きさを表すパラメータであり、Ｗはカーネル関数を表す。）
   請求項１記載の数値計算プログラム。
3. 前記連続関数が、
   

   

   （ｓ_min,a＝ｍ_a／ρ_aＷ（０，ｈ））
   の多項式で表される請求項２記載の数値計算プログラム。
4. 前記連続関数が、
   

   

   である請求項２記載の数値計算プログラム。
5. 連続体を粒子の集合として表現する場合における各粒子について、当該粒子が前記連続体表面にさらされている度合いに対応する係数に応じ且つ放射冷却に基づく内部エネルギーの時間微分を算出し、
   各粒子について算出された前記内部エネルギーの時間微分に基づき、単位時間後の前記内部エネルギーを算出する
   処理を含み、
   前記係数が、
   粒子が前記連続体内部に存在する場合には放射冷却が起こらず、
   粒子が平面上に一様分布している場合には当該粒子と等価な体積を有する、各辺が等長である平行六面体が有する六面のうち一面からの放射冷却が起こり、且つ
   粒子が孤立している場合には当該粒子と等価な体積の球と同じ放射冷却が起こる
   ように補正するための連続関数である、
   コンピュータにより実行される数値計算方法。
6. 連続体を粒子の集合として表現する場合における各粒子について、当該粒子が前記連続体表面にさらされている度合いに対応する係数に応じ且つ放射冷却に基づく内部エネルギーの時間微分を算出する手段と、
   各粒子について算出された前記内部エネルギーの時間微分に基づき、単位時間後の前記内部エネルギーを算出する手段と、
   を有し、
   前記係数が、
   粒子が前記連続体内部に存在する場合には放射冷却が起こらず、
   粒子が平面上に一様分布している場合には当該粒子と等価な体積を有する、各辺が等長である平行六面体が有する六面のうち一面からの放射冷却が起こり、且つ
   粒子が孤立している場合には当該粒子と等価な体積の球と同じ放射冷却が起こる
   ように補正するための連続関数である、
   情報処理装置。

Images