JP5669589B2

JP5669589B2 - 解析装置

Info

Publication number: JP5669589B2
Application number: JP2011009230A
Authority: JP
Inventors: 良孝大西; 大路市嶋
Original assignee: Sumitomo Heavy Industries Ltd
Current assignee: Sumitomo Heavy Industries Ltd
Priority date: 2011-01-19
Filing date: 2011-01-19
Publication date: 2015-02-12
Anticipated expiration: 2031-01-19
Also published as: JP2012150673A

Description

本発明は、粒子系を解析する解析装置およびシミュレーション方法に関する。

従来、古典力学や量子力学等を基に計算機を用いて物質科学全般の現象を探るための方法として、分子動力学法（Molecular Dynamics Method、以下ＭＤ法と称す）に基づくシミュレーションが知られている。ＭＤ法に基づくシミュレーションでは、粒子の運動を物理的により厳密に取り扱えるため相転移等も好適に扱える反面、粒子数が増えると計算量も飛躍的に増大するので実用的にはより少ない数の粒子しか扱えない。したがって従来ＭＤ法は主に材料の物性の予測等、解析対象の形状があまり関係しない用途に使用されることが多かった。

最近、このＭＤ法をマクロスケールの系を扱えるように発展させた繰り込み群分子動力学法（Renormalized Molecular Dynamics、以下ＲＭＤ法と称す）が提案された（例えば、特許文献１、特許文献２参照）。ＲＭＤ法により解析対象は歯車やモータなどの機械的な構造物にまで広がってきている。

解析対象の熱伝導現象について、ＭＤ法やＲＭＤ法では通常、格子振動（フォノン）による熱伝導しか扱えない。したがって自由電子の寄与が大きい金属などの熱伝導については、ＭＤ法やＲＭＤ法では実際とはかけ離れた解析結果が出ることが多い。

これに対して非特許文献１には、空間をブロックで分割し、隣接ブロック同士の温度勾配から、熱伝導を補正する手法が提案されている。

特開２００６−２８５８６６号公報特開２０１０−１７０３０９号公報

大村悦二、福本一郎、宮本勇、「金属を対象とした改良型分子動力学法の開発」、溶接学会論文集、１９９９、第１７巻、第３号、ｐ．３８２−３８６

しかしながら、非特許文献１で提案されている手法では、分割した空間に十分な数の原子が存在しなければ統計平均を取ることができない。また、空間を分割する以上、分割のやり方などの点で複雑な形状を有する解析対象への適用も難しい。

本発明はこうした課題に鑑みてなされたものであり、その目的は、ＭＤ法やＲＭＤ法において熱伝導現象をより正確に扱えるようにする技術の提供にある。

本発明のある態様は解析装置に関する。この解析装置は、複数の粒子を含む系について粒子の運動方程式を数値的に演算することによりその系を解析する解析装置であって、粒子間の距離に基づき粒子に働く力を演算する力演算部と、力演算部によって演算された力と、隣接する粒子間に所与の熱伝導率に基づく熱の流れがある前提で導出された粒子のエネルギの変化量とに基づいて、粒子の位置を演算する粒子位置演算部と、粒子の位置を粒子位置演算部によって演算された位置で更新する位置更新部と、を備える。

この態様によると、隣接する粒子間に所与の熱伝導率に基づく熱の流れがある前提で粒子の位置を演算し、更新することができる。

本発明の別の態様は、シミュレーション方法である。この方法は、分子動力学法を使用して系のふるまいを演算する方法である。この方法では、隣接する粒子間に所与の熱伝導率に基づく熱の流れを仮定し、演算に修正を加える。

なお、以上の構成要素の任意の組み合わせや、本発明の構成要素や表現を装置、方法、システム、コンピュータプログラム、コンピュータプログラムを格納した記録媒体などの間で相互に置換したものもまた、本発明の態様として有効である。

本発明によれば、ＭＤ法やＲＭＤ法において熱伝導現象をより正確に扱うことができる。

実施の形態に係る解析装置の機能および構成を示すブロック図である。図１の解析装置における一連の処理の一例を示すフローチャートである。演算条件設定部によって設定される粒子系の一例を示す模式図である。粒子系の下端からの流入エネルギおよび上端からの流出エネルギを時間の関数として示すグラフである。

以下、本発明を好適な実施の形態をもとに図面を参照しながら説明する。各図面に示される同一または同等の構成要素、部材、処理には、同一の符号を付するものとし、適宜重複した説明は省略する。

ＭＤ法やＲＭＤ法を使用して粒子系を解析する際に、非特許文献１に提案される手法のように粒子系に空間的なアンサンブルを仮定する場合、分割した各空間において意味のある統計平均がとれる程度に十分な数の粒子がその空間に存在することが要求される。しかしながら、ＭＤ法やＲＭＤ法で設定できる系の粒子数は計算に使用する装置の処理能力や計算時間によって制限され、数百から数百万個程度であるのが一般的である。これはアボガドロ数と比べると極めて小さな値である。したがって、分割した各空間に十分な数の粒子を確保することは困難である。

そこで本発明者らは、空間的なアンサンブルの代わりに時間的なアンサンブルを系に導入することに想到した。これにより、空間的なアンサンブルが十分に取れない系、特に複雑な形状の系を解析する場合でも、時間的なアンサンブルを使用することでより正確な物理量を導出することができる。また、そのような正確な物理量が導出される系は、実際の解析対象をより忠実にシミュレート（再現）するものである。

本実施の形態に係る解析装置では、繰り返し演算の各ステップで、系の粒子のもつ運動エネルギを粒子の温度に換算し、その勾配から、フーリエ則にしたがって粒子間でエネルギの受け渡しを行う。各粒子について、このように受け渡されたエネルギ量に基づいて速度が修正される。ここでの熱伝導率（熱伝導係数とも称される）には、実測値や文献値などの陽に定義される値が使用されるので、結果的に自由電子の効果を考慮することができる。これにより、熱に関する物理量をより正確に導出でき、また、特に熱の流れによる寄与が重要な系のシミュレーションを精度良く行うことができる。

図１は、解析装置１００の機能および構成を示すブロック図である。ここに示す各ブロックは、ハードウエア的には、コンピュータのＣＰＵ（central processing unit）をはじめとする素子や機械装置で実現でき、ソフトウエア的にはコンピュータプログラム等によって実現されるが、ここでは、それらの連携によって実現される機能ブロックを描いている。したがって、これらの機能ブロックはハードウエア、ソフトウエアの組合せによっていろいろなかたちで実現できることは、本明細書に触れた当業者には理解されるところである。

解析装置１００は、歯車や梁やモータなどの解析対象を複数の粒子を含む系で記述し、その系について粒子の運動方程式を数値的に演算することによりその系を解析する。特に解析装置１００は演算によって系の時間発展や定常状態を取得し、そのようにして得られたデータから系をシミュレートしたり系の物理量の予測値を提供したりする。

解析装置１００は入力装置１０２および出力装置１０４と接続される。入力装置１０２は、解析装置１００で実行される処理に関係するユーザの入力を受けるためのキーボード、マウスなどであってもよい。入力装置１０２は、インターネットなどのネットワークやＣＤ、ＤＶＤなどの記録媒体から入力を受けるよう構成されていてもよい。出力装置１０４は、ディスプレイなどの表示機器やプリンタなどの印刷機器であってもよい。

解析装置１００は、演算条件設定部１１０と、繰り返し演算部１２０と、結果提示部１５０と、記憶部１６０と、を備える。
演算条件設定部１１０は、解析装置１００における演算に関する演算条件を設定する。演算条件は、解析対象を記述する粒子系Ｓの初期条件と当該解析対象の所与の熱伝導率とを含む。演算条件設定部１１０は、熱伝導率取得部１１２と、モデル生成部１１４と、を含む。

熱伝導率取得部１１２は、入力装置１０２を介してユーザから、解析対象の熱伝導率の情報を取得する。この熱伝導率は公知の文献、例えば「化学便覧」（日本化学会編、丸善株式会社発行）に記載されている文献値であってもよく、あるいはまたユーザの実験による実測値であってもよい。

モデル生成部１１４は、入力装置１０２を介してユーザから取得する入力情報に基づき、Ｎ個の粒子からなる粒子系Ｓを生成する。既に粒子系Ｓが生成され記憶部１６０に記憶されている場合は、モデル生成部１１４はそのように記憶されている情報を記憶部１６０から読み出す。特にモデル生成部１１４は、粒子系ＳのＮ個の粒子の初期位置および初期速度や粒子間の相互作用を表すポテンシャルエネルギ関数や粒子の質量などの以降の演算に必要な初期条件を入力情報から取得し、取得された情報を記憶部１６０に記憶させる。

なお、ミクロスケールの解析対象を記述する場合、ＭＤ法に倣って粒子を原子または分子に対応させてもよい。マクロスケールの解析対象を記述する場合、ＲＭＤ法に倣って粒子を繰り込まれた系の粒子としてもよい。
以下では粒子系Ｓの粒子は全て同質または同等なものとして設定され、かつ、ポテンシャルエネルギ関数は２体のポテンシャルであって粒子によらずに同じ形を有するものとして設定される場合について説明する。しかしながら、他の場合にも本実施の形態に係る技術的思想を適用できることは、本明細書に触れた当業者には明らかである。

繰り返し演算部１２０は、演算条件設定部１１０によって設定された演算条件に基づき所定のアルゴリズムにしたがった繰り返し演算を行う。繰り返し演算部１２０は、力演算部１２２と、粒子位置演算部１２４と、位置更新部１２６と、終了条件判定部１２８と、を含む。以下、繰り返し演算部１２０において繰り返される一連の処理をステップと称す。ひとつのステップは、粒子に働く力の演算と、粒子の位置の演算と、粒子の位置の更新と、終了条件の判定と、を含む。繰り返し演算部１２０においてステップが繰り返されることは、粒子系Ｓの時間が所定の微小な時間刻みΔｔごとに進むことに対応する。

力演算部１２２は記憶部１６０に記憶される粒子系Ｓの情報を参照し、粒子系Ｓの各粒子について、粒子間の距離に基づきその粒子に働く力を演算する。力演算部１２２は、粒子系Ｓのｉ番目（１≦ｉ≦Ｎ）の粒子について、そのｉ番目の粒子との距離が所定のカットオフ距離よりも小さな粒子（以下、近接粒子と称す）を決定する。力演算部１２２は、各近接粒子について、その近接粒子とｉ番目の粒子との間のポテンシャルエネルギ関数およびその近接粒子とｉ番目の粒子との距離に基づいて、その近接粒子がｉ番目の粒子に及ぼす力を演算する。特に力演算部１２２は、その近接粒子とｉ番目の粒子との距離の値におけるポテンシャルエネルギ関数のグラジエント（Gradient）の値から力を算出する。力演算部１２２は、近接粒子がｉ番目の粒子に及ぼす力を全ての近接粒子について足し合わせることによって、ｉ番目の粒子に働く力を算出する。

粒子位置演算部１２４は記憶部１６０に記憶される粒子系Ｓの情報を参照し、粒子系Ｓの各粒子について粒子の位置を演算する。特に粒子位置演算部１２４は、力演算部１２２によって演算された力と、隣接する粒子間に熱伝導率取得部１１２によって取得された所与の熱伝導率に基づく熱の流れがある前提で導出された粒子のエネルギの変化量と、に基づいて、粒子の位置を演算する。粒子位置演算部１２４は、速度算出部１３０と、熱流束密度算出部１３２と、エネルギ変化量算出部１３４と、速度修正部１３６と、位置算出部１３８と、を有する。

速度算出部１３０は、力演算部１２２によって演算された力を含む粒子の運動方程式から粒子の速度を演算する。速度算出部１３０は、粒子系Ｓのｉ番目の粒子について、蛙跳び法やオイラー法などの所定の数値解析の手法に基づき時間刻みΔｔを使用して離散化された粒子の運動方程式に、力演算部１２２によって算出された力を代入することによって、粒子の速度ｖ_ｉ（ベクトル量）を演算する。この演算には以前のステップで演算された粒子の速度が使用される。

熱流束密度算出部１３２は、粒子系Ｓの各粒子に温度が定義される前提で、隣接する粒子間の温度勾配に熱伝導率取得部１１２によって取得された熱伝導率を乗じることによってその粒子間の熱流束密度を演算する。

古典的な気体分子運動論によると、熱平衡状態にある理想気体については以下の式１が成立することが知られている。
ここで、ｍは理想気体の分子（原子の場合もある）の質量、
は理想気体の分子の速度の２乗の平均、ｋはボルツマン定数、Ｔは理想気体の温度である。

本実施の形態では、粒子系Ｓのｉ番目の粒子について、上記式１からの類推で得られる以下の式２を変形し、それによって得られる以下の式３で温度Ｔ_ｉを定義する。
ここでｍ_ｉはｉ番目の粒子の質量である。
なお、式３で定義される粒子の温度を全ての粒子について平均すると粒子系Ｓの温度となるので、式３で定義される粒子の温度を確率変数や温度相当変数と呼んでもよい。

熱流束密度算出部１３２は、粒子系Ｓのｉ番目の粒子について、そのｉ番目の粒子と隣接する粒子（以下、隣接粒子と称す）を決定する。熱流束密度算出部１３２は、各隣接粒子について、その隣接粒子（ｊ番目（ｉ≠ｊ）の粒子とする）とｉ番目の粒子との間の熱流束密度ｑ_ｉｊ（単位は例えばＷ／ｍ^２）をフーリエ則にしたがった以下の式４から算出する。
ここでｒ_ｉｊはｉ番目の粒子とｊ番目の粒子との距離、λは熱伝導率取得部１１２によって取得された熱伝導率である。

エネルギ変化量算出部１３４は、熱流束密度算出部１３２によって演算された熱流束密度に粒子間に定められる面積および速度算出部１３０における演算で使用される時間刻みΔｔを乗じることによって粒子のエネルギの変化量を演算する。

エネルギ変化量算出部１３４は、粒子系Ｓのｉ番目の粒子について、時間刻みΔｔの間に隣接粒子と熱をやりとりすることによってそのｉ番目の粒子から奪われるエネルギΔＥ_ｉ（単位は例えばＪ）を以下の式５から算出する。
ここでΔＳ_ｉｊはｉ番目の粒子とｊ番目の粒子との間に定められる面積であり、例えばＶｏｒｏｎｏｉ解析などを行って求めることができる。ΔＥ_ｉがマイナスの場合はｉ番目の粒子が隣接粒子との熱のやりとりによってエネルギを吸収することを示し、プラスの場合は放出することを示す。

速度修正部１３６は、速度算出部１３０における演算の結果得られた粒子の速度を、エネルギ変化量算出部１３４によって算出された粒子のエネルギの変化量に基づいて変更する。特に速度修正部１３６は、速度算出部１３０における演算の結果得られた粒子の速度によって定まる運動エネルギにその粒子について算出されたエネルギの変化量を加えることで新たな運動エネルギを演算し、演算された新たな運動エネルギを与える粒子の速度を演算する。

本実施の形態では、粒子系Ｓのｉ番目の粒子について、速度算出部１３０における演算の結果得られた粒子の速度ｖ_ｉにより定まる運動エネルギにエネルギ変化量算出部１３４によって算出されたエネルギΔＥ_ｉを加えることによってｉ番目の粒子の運動エネルギを修正する。そのように修正された運動エネルギに対応する速度を新たなｉ番目の粒子の速度とする。すなわち、速度修正部１３６は以下の式６によってｉ番目の粒子の速度を修正する。
ここでｖ'_ｉはｉ番目の粒子の修正された速度である。

位置算出部１３８は、速度修正部１３６によって修正された粒子の速度に基づいて粒子の位置を算出する。位置算出部１３８は、粒子系Ｓのｉ番目の粒子について、所定の数値解析の手法に基づき時間刻みΔｔを使用して離散化された粒子の位置と速度の関係式に、速度修正部１３６によって修正された粒子の速度ｖ'_ｉを適用することによって、粒子の位置ｘ_ｉ（ベクトル量）を演算する。この演算には以前のステップで演算された粒子の位置が使用される。

位置更新部１２６は、記憶部１６０に記憶される粒子系Ｓの各粒子の位置を、位置算出部１３８によって算出された位置で更新する。また位置更新部１２６は、記憶部１６０に記憶される粒子系Ｓの各粒子の速度を、速度修正部１３６によって修正された速度で更新してもよい。

終了条件判定部１２８は、繰り返し演算部１２０における繰り返し演算を終了すべきか否かを判定する。繰り返し演算を終了すべき終了条件は、例えば繰り返し演算が所定の回数行われたことや、粒子系Ｓが定常状態に達したことや、外部から終了の指示を受け付けたことである。終了条件判定部１２８は、終了条件が満たされる場合、繰り返し演算部１２０における繰り返し演算を終了させる。終了条件判定部１２８は、終了条件が満たされない場合、処理を力演算部１２２に戻す。すると力演算部１２２は、位置更新部１２６によって更新された粒子の位置で再び力を演算する。

結果提示部１５０は、演算条件設定部１１０によって設定された粒子系Ｓの解析結果をユーザに提示する。結果提示部１５０は、物理量演算部１５２と、再スケーリング部１５４と、描画制御部１５６と、を含む。

物理量演算部１５２は、繰り返し演算部１２０における繰り返し演算が終了した後、記憶部１６０に記憶される粒子系Ｓの情報に基づき粒子系Ｓの各種物理量、例えば温度や圧力や応力などを演算する。

再スケーリング部１５４は、モデル生成部１１４において粒子系Ｓの粒子が繰り込まれた系の粒子とされた場合、物理量演算部１５２によって演算された物理量を繰り込まれる前の系の物理量に変換する。特に再スケーリング部１５４は、物理量演算部１５２によって演算された物理量に、物理量ごとに定まるスケーリング係数を乗じることによって繰り込まれる前の系の物理量を得る。物理量には応力や温度など繰り込み変換に際して不変となる物理量もあり、そのような物理量についてはスケーリング係数として１（変換前後で不変）が設定される。再スケーリング部１５４は、変換された物理量を出力装置１０４に表示させる。

描画制御部１５６は、繰り返し演算部１２０における繰り返し演算が終了した後、記憶部１６０に記憶される粒子系Ｓの各粒子の位置、速度の情報に基づき、出力装置１０４に粒子系Ｓの時間発展や定常状態の様子をグラフィカルに表示させる。

上述の実施の形態において、記憶部１６０の例は、ハードディスクやメモリである。また、本明細書の記載に基づき、各部を、図示しないＣＰＵや、インストールされたアプリケーションプログラムのモジュールや、システムプログラムのモジュールや、ハードディスクから読み出したデータの内容を一時的に記憶するメモリなどにより実現できることは本明細書に触れた当業者には理解されるところである。

以上の構成による解析装置１００の動作を説明する。
図２は、解析装置１００における一連の処理の一例を示すフローチャートである。熱伝導率取得部１１２は、解析対象に対して設定されるべき熱伝導率を取得する（Ｓ２０２）。モデル生成部１１４は、解析対象を記述した粒子系Ｓを生成する（Ｓ２０４）。力演算部１２２は、粒子間の距離から粒子に働く力を演算する（Ｓ２０６）。速度算出部１３０は、演算された力を含む粒子の運動方程式から速度を演算する（Ｓ２０８）。熱流束密度算出部１３２は、取得された熱伝導率を使用して、粒子間の熱流束密度を演算する（Ｓ２１０）。エネルギ変化量算出部１３４は、演算された熱流束密度に起因する粒子のエネルギの変化量を演算する（Ｓ２１２）。速度修正部１３６は、運動方程式から演算された粒子の速度を、演算された粒子のエネルギの変化量を使用して修正する（Ｓ２１４）。位置算出部１３８は、修正された速度から粒子の位置を算出する（Ｓ２１６）。位置更新部１２６は、記憶部１６０に記憶される粒子の位置を算出された位置で更新する（Ｓ２１８）。終了条件判定部１２８は、終了条件が満たされるか否かを判定する（Ｓ２２０）。終了条件が満たされない場合（Ｓ２２０のＮ）、処理はＳ２０６に戻される。終了条件が満たされる場合（Ｓ２２０のＹ）、結果提示部１５０は演算結果をユーザに提示する（Ｓ２２２）。

本実施の形態に係る解析装置１００では、ＭＤ法やＲＭＤ法を使用して粒子系のふるまいを演算する際、隣接する粒子間に文献や実測等から得られる熱伝導率に基づく熱の流れを仮定し、演算に修正を加える。この熱の流れを演算する際の前提として、各粒子に温度が定義される。この手法では瞬間的に熱エネルギの流れが無秩序となる可能性があるが、解析装置１００では時間刻みΔｔごとに演算を繰り返すので、その繰り返しの回数を適宜増やすことで時間的なアンサンブルを十分に取り、信頼性の高い統計量を導くことができる。

また本実施の形態に係る解析装置１００によると、解析対象の熱伝導、特に従来のＭＤ法やＲＭＤ法では取り扱いが困難な自由電子による熱伝導を補正しながらシミュレーションを行うことができる。その結果、相転移を伴わない通常の解析対象については熱伝導の側面からより正確にシミュレーションを行うことができる。また、一般にＭＤ法やＲＭＤ法に基づく解析は相転移をより忠実に再現できるのであるが、本実施の形態に係る解析装置１００によるとさらに熱伝導現象の再現性も向上するので、結果として液体金属の冷却凝固の過程や固体金属の液化の過程などの熱伝導がクリティカルとなる一次相転移の過程をより正確にシミュレーションできる。

なお、粒子法において、粒子それぞれに温度をパラメータとして与え、熱伝導方程式を解く手法も考えられるが、この場合は潜熱などが取り扱えないので相転移を正確にシミュレートすることができない。これに対して本実施の形態に係る解析装置１００では熱伝導現象および相転移の両方を好適に取り扱える。

また本実施の形態に係る解析装置１００によると、ＲＭＤ法に倣って繰り込まれた粒子系が設定される場合、ミクロスケール、マクロスケールを問わず任意のスケールのシミュレーションが可能となる。

また本実施の形態に係る解析装置１００では、式３、式４、式５、式６など、物理的な根拠がある関係式にしたがって粒子の速度が修正される。これにより、特に熱に関して物理法則に則ったより忠実なシミュレーションが実現される。

図３は、演算条件設定部１１０によって設定される粒子系１０の一例を示す模式図である。ｘ、ｙ、ｚの領域は、それぞれ２．８６×２．４７６×４．４６（nm）である．図３のように、ｚ方向の上端を１００（Ｋ）、下端を２００（Ｋ）の温度で一定となるようにスケーリングを行う。このときスケーリング前後のエネルギ差から流入、流出エネルギを求める。

図４は、粒子系１０の下端からの流入エネルギ１２および上端からの流出エネルギ１４を時間の関数として示すグラフである。図４の縦軸はエネルギ、横軸は時間を示し、この場合の時間は繰り返し演算の回数に比例する。粒子系１０が定常状態になると、図４に示されるように下端からの流入エネルギの絶対値と上端からの流出エネルギの絶対値とが同じになる。

これらのエネルギが求まると、下記の式７から熱伝導率が求まる。
図３、図４の例では、熱伝導率としてアルミニウムの２３８（Ｗ／ｍＫ）を使用した。図４から求めたエネルギの値を式７に代入し、熱伝導率を求めたところ、約２３８±２０．０（Ｗ／ｍＫ）という結果を得たので、熱伝導率について自己無撞着（self-consistent）であることが確認された。
また、本実施の形態に係る手法を用いない、通常のＭＤ法で熱伝導率を求めた場合、値は約３．３（Ｗ／ｍＫ）となった。このことから、本実施の形態に係る手法は目的とする物質の熱伝導率を正しく再現できることが示された。

以上、実施の形態に係る解析装置１００の構成と動作について説明した。これらの実施の形態は例示であり、その各構成要素や各処理の組み合わせにいろいろな変形例が可能なこと、またそうした変形例も本発明の範囲にあることは当業者に理解されるところである。

数値解析の手法にはＶｅｒｌｅｔ法のように、粒子の位置を演算する際に粒子に働く力から粒子の位置を直接演算し、粒子の速度は陽に計算しなくてもよい手法がある。本実施の形態に係る技術的思想をこのような手法に適用する場合、例えば速度修正部は力に時間の次元を有する量を乗算した量に関して上記と同様の修正を行ってもよい。位置算出部は、そのように修正された量に基づいて粒子の位置を演算してもよい。

１００解析装置、１０２入力装置、１０４出力装置、１１０演算条件設定部、１２０繰り返し演算部、１２２力演算部、１２４粒子位置演算部、１２６位置更新部、１５０結果提示部、１６０記憶部。

Claims

複数の粒子を含む系について粒子の運動方程式を数値的に演算することによりその系を解析する解析装置であって、
粒子間の距離に基づき粒子に働く力を演算する力演算部と、
前記力演算部によって演算された力を含む粒子の運動方程式を演算する第１演算部と、前記第１演算部における演算の結果得られる粒子の速度を、隣接する粒子間に所与の熱伝導率に基づく熱の流れがある前提で導出された粒子の運動エネルギから定義された確率変数を仮の温度と定義した熱伝導方程式を解くことによって得られたエネルギの変化量に基づいて変更する第２演算部と、前記第２演算部における変更の結果得られる粒子の速度に基づいて粒子の位置を演算する第３演算部と、を含む粒子位置演算部と、
粒子の位置を前記粒子位置演算部によって演算された位置で更新する位置更新部と、を備えることを特徴とする解析装置。
前記第２演算部は、前記第１演算部における演算の結果得られる粒子の速度によって定まる運動エネルギにその粒子について導出されたエネルギの変化量を加えることで新たな運動エネルギを演算し、演算された新たな運動エネルギを与える粒子の速度を演算し、
前記第３演算部は、前記第２演算部によって演算された粒子の速度に基づいて粒子の位置を演算することを特徴とする請求項１に記載の解析装置。
前記第２演算部は、
各粒子に温度が定義される前提で、隣接する粒子間の温度勾配に所与の熱伝導率を乗じることによってその粒子間の熱流束密度を演算する第４演算部と、
前記第４演算部によって演算された熱流束密度に粒子間に定められる面積および前記第１演算部における演算で使用される時間刻みを乗じることによって粒子のエネルギの変化量を演算する第５演算部と、を有することを特徴とする請求項１または２に記載の解析装置。
複数の粒子を含む系について粒子の運動方程式を数値的に演算することによりその系を解析する機能をコンピュータに実現させるコンピュータプログラムであって、
粒子間の距離に基づき粒子に働く力を演算する機能と、
前記粒子に働く力を演算する機能によって演算された力を含む粒子の運動方程式を演算する機能と、前記粒子の運動方程式を演算する機能による演算の結果得られる粒子の速度を、隣接する粒子間に所与の熱伝導率に基づく熱の流れがある前提で導出された粒子の運動エネルギから定義された確率変数を仮の温度と定義した熱伝導方程式を解くことによって得られたエネルギの変化量に基づいて変更する機能と、前記変更する機能による変更の結果得られる粒子の速度に基づいて粒子の位置を演算する機能と、
粒子の位置を、前記位置を演算する機能により演算された位置で更新する機能と、を前記コンピュータに実現させることを特徴とするコンピュータプログラム。