(57)【要約】【課題】粒子の状態を考慮して放射冷却についての数値計算を可能にする。【解決手段】本方法は、（Ａ）連続体を粒子の集合として表現する場合における各粒子について、当該粒子が連続体表面にさらされている度合いに対応する係数に応じ且つ放射冷却に基づく内部エネルギーの時間微分を算出し、（Ｂ）各粒子について算出された内部エネルギーの時間微分に基づき、単位時間後の内部エネルギーを算出する処理を含む。【選択図】図３
2014211798.tif 000002
【特許請求の範囲】
【請求項１】
連続体を粒子の集合として表現する場合における各粒子について、当該粒子が前記連続体表面にさらされている度合いに対応する係数に応じ且つ放射冷却に基づく内部エネルギーの時間微分を算出し、
各粒子について算出された前記内部エネルギーの時間微分に基づき、単位時間後の前記内部エネルギーを算出する
処理を、コンピュータに実行させるための数値計算プログラム。
【請求項２】
前記係数が、
粒子が前記連続体内部に存在する場合には放射冷却が起こらず、
粒子が平面上に一様分布している場合には当該粒子と等価な体積を有する、各辺が等長である平行六面体が有する六面のうち一面からの放射冷却が起こり、且つ
粒子が孤立している場合には当該粒子と等価な体積の球と同じ放射冷却が起こる
ように補正するための連続関数である
請求項１記載の数値計算プログラム。
【請求項３】
前記連続関数が、規格化された数密度ｓaの関数となっている
【数１】
2014211798.tif 000025
（ｂは、粒子ａの所定範囲内の各粒子を表し、ｍは粒子の質量を表し、ρは密度を表し、ｒabは粒子ａと粒子ｂの距離を表し、ｈは粒子の影響範囲の大きさを表すパラメータであり、Ｗはカーネル関数を表す。）
請求項２記載の数値計算プログラム。
【請求項４】
前記連続関数が、
【数２】
2014211798.tif 000026
（ｓmin,a＝ｍa／ρaＷ（０，ｈ））
の多項式で表される請求項３記載の数値計算プログラム。
【請求項５】
前記連続関数が、
【数３】
2014211798.tif 000027
である請求項３記載の数値計算プログラム。
【請求項６】
連続体を粒子の集合として表現する場合における各粒子について、当該粒子が前記連続体表面にさらされている度合いに対応する係数に応じ且つ放射冷却に基づく内部エネルギーの時間微分を算出し、
各粒子について算出された前記内部エネルギーの時間微分に基づき、単位時間後の前記内部エネルギーを算出する
処理を含み、コンピュータにより実行される数値計算方法。
【請求項７】
連続体を粒子の集合として表現する場合における各粒子について、当該粒子が前記連続体表面にさらされている度合いに対応する係数に応じ且つ放射冷却に基づく内部エネルギーの時間微分を算出する手段と、
各粒子について算出された前記内部エネルギーの時間微分に基づき、単位時間後の前記内部エネルギーを算出する手段と、
を有する情報処理装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、数値計算プログラム、数値計算方法及び情報処理装置に関する。
【背景技術】
【０００２】
数値計算を用いて水や空気等の流体の流れを調べる流体解析や、圧縮されたゴム等の弾性体の振る舞いを調べる弾性体解析において、図１に模式的に示すように、連続体を粒子の分布で表現する手法（以下では総称して粒子法と呼ぶ）が提案されている。そして、粒子法の代表的なものとして、Smoothed Particles Hydrodynamics（ＳＰＨ）法や、Moving Particles Semi-implicit（ＭＰＳ）法が知られている。粒子法における標準的な手法としては、図２に示すように、ある粒子ａに対して予め領域（以下では、「影響領域」と呼ぶ。図２の例では半径２ｈ（ｈは粒子の影響範囲の大きさを表すパラメータ（以下、影響半径と呼ぶ））で規定される円。）を設定しておき、その影響領域の内に存在する相手粒子ｂのみからの相互作用を計算することで物体の運動を解析する。
【０００３】
例えば、鋳造、鍛造などの金属加工では、液体金属の中に冷えて固まった金属が混じる固化過程において体積変化を伴うなど、複雑な過程となる。より具体的には、図３に示すように、容器内部に液体金属を充填すると、容器内部の液体金属の熱は、熱伝導で外側方向に伝わると共に、容器表面及び上面の液体金属表面では放射冷却が発生する。このような現象が進むと、液体金属が固化することになる。図３右側の例では、固化による体積縮小で一部に空間が生じてしまう。
【０００４】
自由表面の扱いの簡単さや複数の計算ノードで並列計算を行う場合の並列性能面、固体との連成計算が比較的容易であるといった利点により鋳造や鍛造シミュレーションは粒子法の活躍が期待される分野であると考えられる。
【０００５】
上で述べたＭＰＳ法については、樹脂流動解析に用いた例が存在している。
【０００６】
また、上で述べたＳＰＨ法を用いた例も存在しており、この例では、複数の粒子の物理量をカーネル関数と呼ぶ重み関数を用いて平滑化して、基礎方程式を離散化している。例えば、鋳造における溶融金属の流れを扱うために、連続を表す式、運動方程式、及びエネルギー方程式を、以下のように離散化している。
【数１】
2014211798.tif 000003
ここで、下付け添え字は、各粒子のインデックスを表す。ρaは粒子ａの密度を表し、ｍaは粒子ａの質量を表し、ｖaは粒子ａの速度ベクトルを表し、ｒaは粒子ａの位置ベクトルを表し、Ｐaは粒子ａの圧力を表し、μaは粒子ａの粘性係数を表し、ｕaは粒子ａの単位質量あたりの内部エネルギーを表し、Ｔaは粒子ａの温度を表し、κaは粒子ａの熱伝導係数を表す。
【０００７】
さらに、以下の関係が成り立つ。（図４参照のこと）
ｖab＝ｖa−ｖb
ｒab（太字）＝ｒa−ｒb
ｒab（イタリック）＝｜ｒab（太字）｜
Ｔab＝Ｔa−Ｔb
【０００８】
また、ηは分母の発散を抑える数値パラメータを表し、ｇは重力加速度を表し、ξは粘性の効果を示す数値パラメータを表し、ξ＝4.96333を用いると良いとされている。Ｗはカーネル関数で、粒子の分布から連続場を構成するのに用いる。例えば、Ｗは、３次のスプライン関数が用いられる。
【０００９】
また、圧力については、以下のような状態方程式を用いて求める。
【数２】
2014211798.tif 000004
但し、Ｐ0はその流体の基準圧力を表し、ρ0は基準密度を表し、γは断熱指数を表す。γ＝１やγ＝７が良く用いられる。
【００１０】
エネルギー方程式である（３）式の右辺は、熱伝導に由来する項である。断熱膨張、圧縮及び粘性による散逸で運動エネルギーが内部エネルギーに転換される効果を考慮する場合には、以下の式を（３）式に加える。
【数３】
2014211798.tif 000005
【００１１】
鋳造などの過程では物体の温度変化を考える上で、図３のような例では、熱伝導と、断熱膨張、圧縮、粘性散逸の項に加えて物体表面からの熱放射による温度変化を考慮することになる。具体的には、黒体の表面から単位面積、単位時間あたりに放出される電磁波のエネルギーが、その黒体の熱力学温度Ｔの４乗に比例するという物理法則であるシュテファン・ボルツマンの法則を考慮に入れることになる。
【００１２】
しかしながら、放射冷却について、粒子の様々な状態を考慮した形で内部エネルギーを算出するような技術は存在していない。具体的には、例えば図５のような状況では、とがった形状の先端（丸印の部分）では、図６のように物体表面が平面となっている場合における当該平面上の粒子（表面粒子と呼ぶ。なお、表面より内側の粒子を非表面粒子又は内部粒子と呼ぶ）よりも放射量が多いはずであり、平面上に存在する粒子よりも単位時間あたりに失われるエネルギーは大きいはずである。しかし、このような違いを数値計算に反映させることが考慮されていない。特に、液滴が飛び散るなどして、孤立粒子になった場合（例えば図７）は全空間に向けて熱放射を行うため、放射冷却が効果的に起きるが、従来の数値計算においてはそのような効果が考慮されていない。
【先行技術文献】
【非特許文献】
【００１３】
【非特許文献１】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）の関係を表す図である。
【図１０】図１０は、本実施の形態に係る情報処理装置の機能ブロック図である。
【図１１】図１１は、本実施の形態に係る処理の処理フローを示す図である。
【図１２】図１２は、本実施の形態に係る処理の処理フローを示す図である。
【図１３】図１３は、コンピュータの機能ブロック図である。
【発明を実施するための形態】
【００１９】
最初に、本発明の実施の形態に係る考え方を述べる。
【００２０】
具体的には、全粒子について、以下の式に従って内部エネルギーの時間微分を算出する。
【数４】
2014211798.tif 000006
【００２１】
ここでｇaは、粒子ａが物体表面（連続体表面）にさらされている度合いを定量化するための係数であって、
（Ａ）粒子が連続体内部に存在する場合には放射冷却が起こらず、
（Ｂ）粒子が平面上に一様分布している場合には放射冷却が、当該粒子と等価な体積を有する、各辺が等長である平行六面体が有する六面のうち一面からの放射冷却が起こり、且つ
（Ｃ）粒子が孤立している場合には当該粒子と等価な体積の球と同じ放射冷却が起こる
ように補正するための連続関数である。なお、シュテファン・ボルツマンの法則は、（８）式において、係数ｇa以外の部分で表される。また、σは、物質に依存しないシュテファン・ボルツマン係数であり、αは物体の、黒体からのずれを表す係数（以下、放射率と呼ぶ）を表し、ｄは空間次元を表す。ｄは、１、２又は３である。さらに、（ｍa／ρa）1-(1/d)は、粒子ａの表面積を表している。ここで、「当該粒子と等価な体積を有する、各辺が等長である平行六面体が有する六面のうち一面」は、当該粒子と等価な体積の立方体が有する六面のうち一面と等しい面積を有することになる。また、係数ｇaの満たすべき（Ｂ）の条件の代わりとして、表面形状に応じて粒子を、当該粒子と等価な体積の多面体や半球と仮定した場合において、当該粒子と等価な体積の立方体が有する六面のうち一面と等しい面積を有する立体であっても良い。
【００２２】
本実施の形態では、このため、係数ｇaを、以下に示すような、規格化された数密度の関数として表す。
【数５】
2014211798.tif 000007
【００２３】
この規格化された数密度は、着目している粒子ａの影響範囲内に他の粒子が十分にあると１に近い値となり、影響範囲内の粒子が少ない場合には小さな値となる。
【００２４】
ここで、カーネル関数Ｗ（ｒ，ｈ）は、例えば以下のような５次のスプライン関数である。
【数６】
2014211798.tif 000008
ここでαdは規格化係数である。また、ｑ＝ｒ／ｈである。
【００２５】
規格化された数密度ｓaは、物体内部では「１」、孤立した粒子では以下の値となる。
【数７】
2014211798.tif 000009
【００２６】
さらに、平面とみなせる平面上では「（１＋ｓmin,a）／２」という値になるように定義される。
【００２７】
このような条件を満たすような係数ｇaは一意には決まらないが、例えば、以下のような多項式を採用すれば、上で述べた（Ａ）乃至（Ｃ）の条件は満たされる。
【数８】
2014211798.tif 000010
ここで、記号は以下のような意味である。
【数９】
2014211798.tif 000011
【００２８】
以下、この式で（Ａ）乃至（Ｃ）の条件を満たしていることを説明する。
【００２９】
Ｑ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を除外した式の値）に一致するためである。
（Ｃ）孤立した粒子は等価な体積の球と同様に放射冷却が生ずる。
この条件は、以下のような式で表される。
【数１０】
2014211798.tif 000012
ここでＳ1,aは、上で述べたとおりに定義され、粒子ａを球とみなしたときの表面積を表す。（ｍa／ρa）1-(1/d)は平面上に一様に分布している配置での粒子ａ１つあたりの表面積を表す。孤立した粒子については、これに係数ｇaを掛けた結果が、粒子ａを球とみなした時の表面積と一致させたならば、等価な体積の球と同様に放射冷却をするという条件を満たすことになる。
【００３１】
ここで、数値計算上、内部粒子についてはＱaの値が粒子配置の多少のずれによって０付近で多少揺らぐと考えられるが、それによって放射冷却の補正係数が敏感に変動することを避けることが好ましい。従って、下記の式の条件を付加する。
【数１１】
2014211798.tif 000013
【００３２】
よって、以下に示す４式を満たし且つ最も次数の低い多項式を求めると（９）式となる。
【数１２】
2014211798.tif 000014
【００３３】
但し、最も次数の低い多項式ではなく、図９に示すような滑らかな関数となれば、より次数の高い多項式などを採用するようにしても良い。
【００３４】
なお、３次元においてはＡa＝３2/3（４π）1/3となるので、この時のＱaとｇa（Ｑa）との関係は、図９に示すようになる。
【００３５】
図９では、Ｑa＝０付近ではｇa（Ｑa）は滑らかに増加し、Ｑa＝１／２ではｇa（Ｑa）＝１となっている。Ｑaがそれより大きくなると、ｇa（Ｑa）は急激に増加する。
【００３６】
このような係数ｇaを採用することによって、粒子の状態を放射冷却による内部エネルギーの減少に適切に反映させることができる。また、このような係数ｇaは、規格化された数密度を変数とする多項式で表される連続関数であって、処理フローにおける分岐といった形ではなく無段階で粒子の状態を反映させることができるようになっている。
【００３７】
次に、このような放射冷却に係るシミュレーションを行うための情報処理装置１００の機能ブロック図を図１０に示す。
【００３８】
本実施の形態に係る情報処理装置１００は、入力部１１０と、第１データ格納部１２０と、物理量算出部１３０と、第２データ格納部１４０と、出力部１５０とを有する。
【００３９】
入力部１１０は、例えばネットワークで接続されている他のコンピュータからデータを取得して又はユーザからのデータ入力を受け付けて、処理すべきデータとして第１データ格納部１２０に格納する。
【００４０】
入力データは、数値計算を行う対象である連続体粒子のデータと、連続体粒子の運動に関する境界条件などを設定する固定境界要素のデータとを含む。連続体粒子は、流体をモデル化したものであって、これを表現するためのデータは、例えば初期の中心位置座標、初期の速度、影響半径、密度、質量、粘性等である。固定境界要素は、鋳型などの表面などを微小部分に分割した面要素をモデル化したものであって、このデータとしては、例えば微小な円盤の集合として境界全体をモデル化して、各境界要素に中心座標、面要素の法線ベクトル、面要素の面積を含む。ポリゴンの集合として境界全体を表現して、各境界要素に複数の頂点の位置座標を設定するようにしても良い。
【００４１】
物理量算出部１３０は、近傍リスト生成部１３１と、放射冷却算出部１３２と、積分処理部１３３とを有する。物理量算出部１３０は、単位時間毎に各粒子の物理量を算出する処理を行う。そのために、近傍リスト生成部１３１は、単位時間毎に、各粒子の影響範囲に入っている他の粒子のリストを生成し、第２データ格納部１４０に格納する。また、放射冷却算出部１３２は、単位時間毎に、各粒子について放射冷却による内部エネルギーの時間微分を算出し、第２データ格納部１４０に格納する。積分処理部１３３は、各粒子について算出された加速度、速度、密度の時間微分、内部エネルギーの時間微分等を時間積分して、１単位時間後の物理量（例えば速度、位置、密度、内部エネルギー等）を算出し、第２データ格納部１４０に格納する。
【００４２】
出力部１５０は、第２データ格納部１４０に格納されている単位時間毎の物理量を用いて出力データを生成し、他のコンピュータに出力したり、印刷装置や表示装置などの出力装置に出力する。
【００４３】
次に、図１１及び図１２を用いて情報処理装置１００の処理内容について説明する。まず、入力部１１０は、他のコンピュータから処理に用いられるデータを取得するか、ユーザからの入力を受け付けて、第１データ格納部１２０に格納する（ステップＳ１）。上で述べたような入力データを取得して、第１データ格納部１２０に格納する。
【００４４】
次に、物理量算出部１３０は、時刻ｔを０に初期化し（ステップＳ３）、近傍リスト生成部１３１は、時刻ｔにおける粒子分布に基づいて、各粒子について近傍粒子リストを生成し、第２データ格納部１４０に格納する（ステップＳ５）。
【００４５】
ｔ＝０の場合には、第１データ格納部１２０に格納されているデータに含まれる初期位置から、例えば着目している粒子から距離が影響範囲２ｈ以内の他の粒子の識別子を、リストに追加する。影響範囲は、相互に影響を与え合う範囲であり、例えば、粒子が移動する場合に、他の粒子へ力を加える等の処理を行う範囲である。なお、ｔ＝１以降については、第２データ格納部１４０に格納されているデータを用いて、近傍粒子リストを生成する。
【００４６】
そして、物理量算出部１３０は、未処理の粒子ａを１つ特定する（ステップＳ７）。その後、物理量算出部１３０は、連続体（物体）の物理モデル（例えば（１）式乃至（３）式）に応じて、特定された粒子ａについての解析対象量（放射冷却以外の要素に由来する内部エネルギーの時間微分を含む）を、当該粒子及びその近傍粒子リストに含まれる粒子との相互作用の重ね合わせとして算出し、第２データ格納部１４０に格納する（ステップＳ９）。このステップの処理については、例えば非特許文献２などで説明されており、ここでは説明を省略する。
【００４７】
解析対象量としては、例えば加速度、速度、密度の時間微分、放射冷却以外の要素に由来する内部エネルギーの時間微分が含まれる。放射冷却以外の要素に由来する内部エネルギーの時間微分は、以下のように表される。
【数１３】
2014211798.tif 000015
ｕaは、粒子ａの内部エネルギーを表す。
【００４８】
さらに、放射冷却算出部１３２は、特定された粒子ａについて、物体表面にさらされている度合いに応じた補正係数ｇaを用いて、放射冷却によって単位時間あたりに失われる内部エネルギー量（内部エネルギーの時間微分）を算出し、第２データ格納部１４０に格納する（ステップＳ１１）。具体的には、（８）式及び（９）式に従って、内部エネルギーの時間微分を算出する。ここでは、放射冷却に由来する内部エネルギーの時間微分は、以下のように表される。
【数１４】
2014211798.tif 000016
【００４９】
そして、物理量算出部１３０は、算出された内部エネルギー量を、内部エネルギーの時間微分に加算し、第２データ格納部１４０に格納する（ステップＳ１３）。具体的には、以下のような演算を行う。
【数１５】
2014211798.tif 000017
【００５０】
そして、物理量算出部１３０は、未処理の粒子が存在するか判断する（ステップＳ１５）。未処理の粒子が存在する場合には、処理はステップＳ７に戻る。一方、未処理の粒子が存在しない場合には、処理は端子Ａを介して図１２の処理に移行する。
【００５１】
図１２の処理の説明に移行して、物理量算出部１３０の積分処理部１３３は、各粒子について、解析対象量の時間積分を実行し、処理結果を第２データ格納部１４０に格納する（ステップＳ１７）。着目している内部エネルギーについては、以下のような演算を実行する。
【数１６】
2014211798.tif 000018
ここでｕtaは時刻ｔにおける内部エネルギーを表している。
【００５２】
その他の解析対象量についても、同様の演算を行って、各粒子について、時刻ｔ＋１についての物理量を算出する。
【００５３】
そして、出力部１５０は、第２データ格納部１４０に格納されている時刻ｔ＋１についての物理量を、出力装置（他のコンピュータ、印刷装置又は表示装置等）に出力する（ステップＳ１９）。
【００５４】
物理量算出部１３０は、例えば時刻ｔが処理終了時刻になったか否かを判断することで処理終了であるか否かを判断する（ステップＳ２１）。処理終了であれば、処理を終了し、処理終了でなければ、物理量算出部１３０は、時刻ｔを１インクリメントして（ステップＳ２３）、処理は端子Ｂを介してステップＳ５に移行する。
【００５５】
以上のような処理を行えば、係数ｇaは、粒子毎にその状態を反映した値、すなわち、表面にさらされている度合いに応じた値となる。そして、各粒子が、単位時間毎に異なる状態であっても、その状態に応じた放射冷却による内部エネルギーの時間微分が算出できるようになる。
【００５６】
この際、処理フローの分岐という形ではなく、連続関数として補正係数ｇaが定義されているので、放射冷却の効果を無段階に反映させることができ、正確に物体の温度変化を計算することができるようになる。
【００５７】
以上本発明の実施の形態を説明したが、本発明はこれに限定されるものではない。例えば、図１０に示した機能ブロック図は一例であって、プログラムモジュール構成とは一致しない場合がある。また、処理フローについても、各粒子の演算については複数のプロセッサで並列演算を行うように処理フローを変形したりするようにしても良い。
【００５８】
また、情報処理装置１００は、１台のコンピュータではなく、複数台のコンピュータでその機能が分担される場合もある。
【００５９】
なお、上で述べた情報処理装置１００は、コンピュータ装置であって、図１３に示すように、メモリ２５０１とＣＰＵ（Central Processing Unit）２５０３とハードディスク・ドライブ（ＨＤＤ：Hard Disk Drive）２５０５と表示装置２５０９に接続される表示制御部２５０７とリムーバブル・ディスク２５１１用のドライブ装置２５１３と入力装置２５１５とネットワークに接続するための通信制御部２５１７とがバス２５１９で接続されている。オペレーティング・システム（ＯＳ：Operating System）及び本実施例における処理を実施するためのアプリケーション・プログラムは、ＨＤＤ２５０５に格納されており、ＣＰＵ２５０３により実行される際にはＨＤＤ２５０５からメモリ２５０１に読み出される。ＣＰＵ２５０３は、アプリケーション・プログラムの処理内容に応じて表示制御部２５０７、通信制御部２５１７、ドライブ装置２５１３を制御して、所定の動作を行わせる。また、処理途中のデータについては、主としてメモリ２５０１に格納されるが、ＨＤＤ２５０５に格納されるようにしてもよい。本技術の実施例では、上で述べた処理を実施するためのアプリケーション・プログラムはコンピュータ読み取り可能なリムーバブル・ディスク２５１１に格納されて頒布され、ドライブ装置２５１３からＨＤＤ２５０５にインストールされる。インターネットなどのネットワーク及び通信制御部２５１７を経由して、ＨＤＤ２５０５にインストールされる場合もある。このようなコンピュータ装置は、上で述べたＣＰＵ２５０３、メモリ２５０１などのハードウエアとＯＳ及びアプリケーション・プログラムなどのプログラムとが協働することにより、上で述べたような各種機能を実現する。
【００６０】
以上述べた本実施の形態をまとめると、以下のようになる。
【００６１】
本実施の形態に係る数値計算方法は、（Ａ）連続体を粒子の集合として表現する場合における各粒子について、当該粒子が連続体表面にさらされている度合いに対応する係数に応じ且つ放射冷却に基づく内部エネルギーの時間微分を算出し、（Ｂ）各粒子について算出された内部エネルギーの時間微分に基づき、単位時間後の内部エネルギーを算出する処理を含む。
【００６２】
このようにすれば、個々の粒子の状態を反映させることができ、状態変化が生じてもその状態変化に応じて、放射冷却に基づき内部エネルギーの時間微分を算出できるようになる。
【００６３】
上で述べた係数が、（ａ）粒子が連続体内部に存在する場合には放射冷却が起こらず、（ｂ）粒子が平面上に一様分布している場合には当該粒子と等価な体積の立方体のうち一面からの放射冷却が起こり、且つ（ｃ）粒子が孤立している場合には当該粒子と等価な体積の球と同じ放射冷却が起こるように補正するための連続関数である場合もある。
【００６４】
このようにすれば放射冷却による内部エネルギーの時間微分を正確に算出できるようになる。
【００６５】
さらに、上で述べた連続関数が、規格化された数密度ｓaの関数となっている
【数１７】
2014211798.tif 000019
（ｂは、粒子ａの所定範囲内の各粒子を表し、ｍは粒子の質量を表し、ρは密度を表し、ｒabは粒子ａと粒子ｂの距離を表し、ｈは粒子の影響範囲の大きさを表すパラメータであり、Ｗはカーネル関数を表す。）場合もある。規格化された数密度は、影響範囲内の粒子ｂに応じて変化するため、粒子ａの状態に応じた放射冷却による内部エネルギーの時間微分が算出されるようになる。
【００６６】
さらに、上で述べた連続関数が、
【数１８】
2014211798.tif 000020
（ｓmin,a＝ｍa／ρaＷ（０，ｈ））
の多項式で表される場合もある。このような変数変換を行うことで、多項式についての条件設定が容易になる。
【００６７】
さらに、上で述べた連続関数が、
【数１９】
2014211798.tif 000021
である場合もある。上で述べた（ａ）乃至（ｃ）の３条件に対応する条件とＱa＝０付近で滑らかであるという条件とを満たす次数の最も小さいものを採用することで、滑らかに変化する連続関数を用いることができるようになる。
【００６８】
なお、上で述べたような処理をコンピュータに実行させるためのプログラムを作成することができ、当該プログラムは、例えばフレキシブル・ディスク、ＣＤ−ＲＯＭなどの光ディスク、光磁気ディスク、半導体メモリ（例えばＲＯＭ）、ハードディスク等のコンピュータ読み取り可能な記憶媒体又は記憶装置に格納される。なお、処理途中のデータについては、ＲＡＭ等の記憶装置に一時保管される。
【００６９】
以上の実施例を含む実施形態に関し、さらに以下の付記を開示する。
【００７０】
（付記１）
連続体を粒子の集合として表現する場合における各粒子について、当該粒子が前記連続体表面にさらされている度合いに対応する係数に応じ且つ放射冷却に基づく内部エネルギーの時間微分を算出し、
各粒子について算出された前記内部エネルギーの時間微分に基づき、単位時間後の前記内部エネルギーを算出する
処理を、コンピュータに実行させるための数値計算プログラム。
【００７１】
（付記２）
前記係数が、
粒子が前記連続体内部に存在する場合には放射冷却が起こらず、
粒子が平面上に一様分布している場合には当該粒子と等価な体積を有する、各辺が等長である平行六面体が有する六面のうち一面からの放射冷却が起こり、且つ
粒子が孤立している場合には当該粒子と等価な体積の球と同じ放射冷却が起こる
ように補正するための連続関数である
付記１記載の数値計算プログラム。
【００７２】
（付記３）
前記連続関数が、規格化された数密度ｓaの関数となっている
【数２０】
2014211798.tif 000022
（ｂは、粒子ａの所定範囲内の各粒子を表し、ｍは粒子の質量を表し、ρは密度を表し、ｒabは粒子ａと粒子ｂの距離を表し、ｈは粒子の影響範囲の大きさを表すパラメータであり、Ｗはカーネル関数を表す。）
付記２記載の数値計算プログラム。
【００７３】
（付記４）
前記連続関数が、
【数２１】
2014211798.tif 000023
（ｓmin,a＝ｍa／ρaＷ（０，ｈ））
の多項式で表される付記３記載の数値計算プログラム。
【００７４】
（付記５）
前記連続関数が、
【数２２】
2014211798.tif 000024
である付記３記載の数値計算プログラム。
【００７５】
（付記６）
連続体を粒子の集合として表現する場合における各粒子について、当該粒子が前記連続体表面にさらされている度合いに対応する係数に応じ且つ放射冷却に基づく内部エネルギーの時間微分を算出し、
各粒子について算出された前記内部エネルギーの時間微分に基づき、単位時間後の前記内部エネルギーを算出する
処理を含み、コンピュータにより実行される数値計算方法。
【００７６】
（付記７）
連続体を粒子の集合として表現する場合における各粒子について、当該粒子が前記連続体表面にさらされている度合いに対応する係数に応じ且つ放射冷却に基づく内部エネルギーの時間微分を算出する手段と、
各粒子について算出された前記内部エネルギーの時間微分に基づき、単位時間後の前記内部エネルギーを算出する手段と、
を有する情報処理装置。
【符号の説明】
【００７７】
１１０  入力部
１２０  第１データ格納部
１３０  物理量算出部
１４０  第２データ格納部
１５０  出力部
１３１  近傍リスト生成部
１３２  放射冷却算出部
１３３  積分処理部
【図１】
2014211798.tif 000028
【図２】
2014211798.tif 000029
【図３】
2014211798.tif 000030
【図４】
2014211798.tif 000031
【図５】
2014211798.tif 000032
【図６】
2014211798.tif 000033
【図７】
2014211798.tif 000034
【図８】
2014211798.tif 000035
【図９】
2014211798.tif 000036
【図１０】
2014211798.tif 000037
【図１１】
2014211798.tif 000038
【図１２】
2014211798.tif 000039
【図１３】
2014211798.tif 000040