複数の粒子を含む系について粒子の運動方程式を数値的に演算することによりその系を解析する解析装置であって、粒子間の距離に基づき粒子に働く力を演算する力演算部と、前記力演算部によって演算された力と、隣接する粒子間に所与の熱伝導率に基づく熱の流れがある前提で導出された粒子のエネルギの変化量とに基づいて、粒子の位置を演算する粒子位置演算部と、粒子の位置を前記粒子位置演算部によって演算された位置で更新する位置更新部と、を備えることを特徴とする解析装置。前記粒子位置演算部は、前記力演算部によって演算された力を含む粒子の運動方程式を演算する第1演算部と、前記第1演算部における演算の結果得られる粒子の速度を、隣接する粒子間に所与の熱伝導率に基づく熱の流れがある前提で導出された粒子のエネルギの変化量に基づいて変更する第2演算部と、前記第2演算部における変更の結果得られる粒子の速度に基づいて粒子の位置を演算する第3演算部と、を含むことを特徴とする請求項1に記載の解析装置。前記第2演算部は、前記第1演算部における演算の結果得られる粒子の速度によって定まる運動エネルギにその粒子について導出されたエネルギの変化量を加えることで新たな運動エネルギを演算し、演算された新たな運動エネルギを与える粒子の速度を演算し、前記第3演算部は、前記第2演算部によって演算された粒子の速度に基づいて粒子の位置を演算することを特徴とする請求項2に記載の解析装置。前記第2演算部は、各粒子に温度が定義される前提で、隣接する粒子間の温度勾配に所与の熱伝導率を乗じることによってその粒子間の熱流束密度を演算する第4演算部と、前記第4演算部によって演算された熱流束密度に粒子間に定められる面積および前記第1演算部における演算で使用される時間刻みを乗じることによって粒子のエネルギの変化量を演算する第5演算部と、を有することを特徴とする請求項2または3に記載の解析装置。分子動力学法を使用して系のふるまいを演算する方法において、隣接する粒子間に所与の熱伝導率に基づく熱の流れを仮定し、演算に修正を加えたことを特徴とするシミュレーション方法。複数の粒子を含む系について粒子の運動方程式を数値的に演算することによりその系を解析する機能をコンピュータに実現させるコンピュータプログラムであって、粒子間の距離に基づき粒子に働く力を演算する機能と、演算された力と、隣接する粒子間に所与の熱伝導率に基づく熱の流れがある前提で導出された粒子のエネルギの変化量とに基づいて、粒子の位置を演算する機能と、粒子の位置を、演算された位置で更新する機能と、を前記コンピュータに実現させることを特徴とするコンピュータプログラム。しかしながら、非特許文献1で提案されている手法では、分割した空間に十分な数の原子が存在しなければ統計平均を取ることができない。また、空間を分割する以上、分割のやり方などの点で複雑な形状を有する解析対象への適用も難しい。本発明はこうした課題に鑑みてなされたものであり、その目的は、MD法やRMD法において熱伝導現象をより正確に扱えるようにする技術の提供にある。本発明のある態様は解析装置に関する。この解析装置は、複数の粒子を含む系について粒子の運動方程式を数値的に演算することによりその系を解析する解析装置であって、粒子間の距離に基づき粒子に働く力を演算する力演算部と、力演算部によって演算された力と、隣接する粒子間に所与の熱伝導率に基づく熱の流れがある前提で導出された粒子のエネルギの変化量とに基づいて、粒子の位置を演算する粒子位置演算部と、粒子の位置を粒子位置演算部によって演算された位置で更新する位置更新部と、を備える。この態様によると、隣接する粒子間に所与の熱伝導率に基づく熱の流れがある前提で粒子の位置を演算し、更新することができる。本発明の別の態様は、シミュレーション方法である。この方法は、分子動力学法を使用して系のふるまいを演算する方法である。この方法では、隣接する粒子間に所与の熱伝導率に基づく熱の流れを仮定し、演算に修正を加える。なお、以上の構成要素の任意の組み合わせや、本発明の構成要素や表現を装置、方法、システム、コンピュータプログラム、コンピュータプログラムを格納した記録媒体などの間で相互に置換したものもまた、本発明の態様として有効である。本発明によれば、MD法やRMD法において熱伝導現象をより正確に扱うことができる。以下、本発明を好適な実施の形態をもとに図面を参照しながら説明する。各図面に示される同一または同等の構成要素、部材、処理には、同一の符号を付するものとし、適宜重複した説明は省略する。MD法やRMD法を使用して粒子系を解析する際に、非特許文献1に提案される手法のように粒子系に空間的なアンサンブルを仮定する場合、分割した各空間において意味のある統計平均がとれる程度に十分な数の粒子がその空間に存在することが要求される。しかしながら、MD法やRMD法で設定できる系の粒子数は計算に使用する装置の処理能力や計算時間によって制限され、数百から数百万個程度であるのが一般的である。これはアボガドロ数と比べると極めて小さな値である。したがって、分割した各空間に十分な数の粒子を確保することは困難である。そこで本発明者らは、空間的なアンサンブルの代わりに時間的なアンサンブルを系に導入することに想到した。これにより、空間的なアンサンブルが十分に取れない系、特に複雑な形状の系を解析する場合でも、時間的なアンサンブルを使用することでより正確な物理量を導出することができる。また、そのような正確な物理量が導出される系は、実際の解析対象をより忠実にシミュレート再現するものである。本実施の形態に係る解析装置では、繰り返し演算の各ステップで、系の粒子のもつ運動エネルギを粒子の温度に換算し、その勾配から、フーリエ則にしたがって粒子間でエネルギの受け渡しを行う。各粒子について、このように受け渡されたエネルギ量に基づいて速度が修正される。ここでの熱伝導率熱伝導係数とも称されるには、実測値や文献値などの陽に定義される値が使用されるので、結果的に自由電子の効果を考慮することができる。これにより、熱に関する物理量をより正確に導出でき、また、特に熱の流れによる寄与が重要な系のシミュレーションを精度良く行うことができる。図1は、解析装置100の機能および構成を示すブロック図である。ここに示す各ブロックは、ハードウエア的には、コンピュータのCPUcentralprocessingunitをはじめとする素子や機械装置で実現でき、ソフトウエア的にはコンピュータプログラム等によって実現されるが、ここでは、それらの連携によって実現される機能ブロックを描いている。したがって、これらの機能ブロックはハードウエア、ソフトウエアの組合せによっていろいろなかたちで実現できることは、本明細書に触れた当業者には理解されるところである。解析装置100は、歯車や梁やモータなどの解析対象を複数の粒子を含む系で記述し、その系について粒子の運動方程式を数値的に演算することによりその系を解析する。特に解析装置100は演算によって系の時間発展や定常状態を取得し、そのようにして得られたデータから系をシミュレートしたり系の物理量の予測値を提供したりする。解析装置100は入力装置102および出力装置104と接続される。入力装置102は、解析装置100で実行される処理に関係するユーザの入力を受けるためのキーボード、マウスなどであってもよい。入力装置102は、インターネットなどのネットワークやCD、DVDなどの記録媒体から入力を受けるよう構成されていてもよい。出力装置104は、ディスプレイなどの表示機器やプリンタなどの印刷機器であってもよい。解析装置100は、演算条件設定部110と、繰り返し演算部120と、結果提示部150と、記憶部160と、を備える。演算条件設定部110は、解析装置100における演算に関する演算条件を設定する。演算条件は、解析対象を記述する粒子系Sの初期条件と当該解析対象の所与の熱伝導率とを含む。演算条件設定部110は、熱伝導率取得部112と、モデル生成部114と、を含む。熱伝導率取得部112は、入力装置102を介してユーザから、解析対象の熱伝導率の情報を取得する。この熱伝導率は公知の文献、例えば化学便覧日本化学会編、丸善株式会社発行に記載されている文献値であってもよく、あるいはまたユーザの実験による実測値であってもよい。モデル生成部114は、入力装置102を介してユーザから取得する入力情報に基づき、N個の粒子からなる粒子系Sを生成する。既に粒子系Sが生成され記憶部160に記憶されている場合は、モデル生成部114はそのように記憶されている情報を記憶部160から読み出す。特にモデル生成部114は、粒子系SのN個の粒子の初期位置および初期速度や粒子間の相互作用を表すポテンシャルエネルギ関数や粒子の質量などの以降の演算に必要な初期条件を入力情報から取得し、取得された情報を記憶部160に記憶させる。なお、ミクロスケールの解析対象を記述する場合、MD法に倣って粒子を原子または分子に対応させてもよい。マクロスケールの解析対象を記述する場合、RMD法に倣って粒子を繰り込まれた系の粒子としてもよい。以下では粒子系Sの粒子は全て同質または同等なものとして設定され、かつ、ポテンシャルエネルギ関数は2体のポテンシャルであって粒子によらずに同じ形を有するものとして設定される場合について説明する。しかしながら、他の場合にも本実施の形態に係る技術的思想を適用できることは、本明細書に触れた当業者には明らかである。繰り返し演算部120は、演算条件設定部110によって設定された演算条件に基づき所定のアルゴリズムにしたがった繰り返し演算を行う。繰り返し演算部120は、力演算部122と、粒子位置演算部124と、位置更新部126と、終了条件判定部128と、を含む。以下、繰り返し演算部120において繰り返される一連の処理をステップと称す。ひとつのステップは、粒子に働く力の演算と、粒子の位置の演算と、粒子の位置の更新と、終了条件の判定と、を含む。繰り返し演算部120においてステップが繰り返されることは、粒子系Sの時間が所定の微小な時間刻みΔtごとに進むことに対応する。力演算部122は記憶部160に記憶される粒子系Sの情報を参照し、粒子系Sの各粒子について、粒子間の距離に基づきその粒子に働く力を演算する。力演算部122は、粒子系Sのi番目1≦i≦Nの粒子について、そのi番目の粒子との距離が所定のカットオフ距離よりも小さな粒子以下、近接粒子と称すを決定する。力演算部122は、各近接粒子について、その近接粒子とi番目の粒子との間のポテンシャルエネルギ関数およびその近接粒子とi番目の粒子との距離に基づいて、その近接粒子がi番目の粒子に及ぼす力を演算する。特に力演算部122は、その近接粒子とi番目の粒子との距離の値におけるポテンシャルエネルギ関数のグラジエントGradientの値から力を算出する。力演算部122は、近接粒子がi番目の粒子に及ぼす力を全ての近接粒子について足し合わせることによって、i番目の粒子に働く力を算出する。粒子位置演算部124は記憶部160に記憶される粒子系Sの情報を参照し、粒子系Sの各粒子について粒子の位置を演算する。特に粒子位置演算部124は、力演算部122によって演算された力と、隣接する粒子間に熱伝導率取得部112によって取得された所与の熱伝導率に基づく熱の流れがある前提で導出された粒子のエネルギの変化量と、に基づいて、粒子の位置を演算する。粒子位置演算部124は、速度算出部130と、熱流束密度算出部132と、エネルギ変化量算出部134と、速度修正部136と、位置算出部138と、を有する。速度算出部130は、力演算部122によって演算された力を含む粒子の運動方程式から粒子の速度を演算する。速度算出部130は、粒子系Sのi番目の粒子について、蛙跳び法やオイラー法などの所定の数値解析の手法に基づき時間刻みΔtを使用して離散化された粒子の運動方程式に、力演算部122によって算出された力を代入することによって、粒子の速度viベクトル量を演算する。この演算には以前のステップで演算された粒子の速度が使用される。熱流束密度算出部132は、粒子系Sの各粒子に温度が定義される前提で、隣接する粒子間の温度勾配に熱伝導率取得部112によって取得された熱伝導率を乗じることによってその粒子間の熱流束密度を演算する。古典的な気体分子運動論によると、熱平衡状態にある理想気体については以下の式1が成立することが知られている。ここで、mは理想気体の分子原子の場合もあるの質量、は理想気体の分子の速度の2乗の平均、kはボルツマン定数、Tは理想気体の温度である。本実施の形態では、粒子系Sのi番目の粒子について、上記式1からの類推で得られる以下の式2を変形し、それによって得られる以下の式3で温度Tiを定義する。ここでmiはi番目の粒子の質量である。なお、式3で定義される粒子の温度を全ての粒子について平均すると粒子系Sの温度となるので、式3で定義される粒子の温度を確率変数や温度相当変数と呼んでもよい。熱流束密度算出部132は、粒子系Sのi番目の粒子について、そのi番目の粒子と隣接する粒子以下、隣接粒子と称すを決定する。熱流束密度算出部132は、各隣接粒子について、その隣接粒子j番目i≠jの粒子とするとi番目の粒子との間の熱流束密度qij単位は例えばW/m2をフーリエ則にしたがった以下の式4から算出する。ここでrijはi番目の粒子とj番目の粒子との距離、λは熱伝導率取得部112によって取得された熱伝導率である。エネルギ変化量算出部134は、熱流束密度算出部132によって演算された熱流束密度に粒子間に定められる面積および速度算出部130における演算で使用される時間刻みΔtを乗じることによって粒子のエネルギの変化量を演算する。エネルギ変化量算出部134は、粒子系Sのi番目の粒子について、時間刻みΔtの間に隣接粒子と熱をやりとりすることによってそのi番目の粒子から奪われるエネルギΔEi単位は例えばJを以下の式5から算出する。ここでΔSijはi番目の粒子とj番目の粒子との間に定められる面積であり、例えばVoronoi解析などを行って求めることができる。ΔEiがマイナスの場合はi番目の粒子が隣接粒子との熱のやりとりによってエネルギを吸収することを示し、プラスの場合は放出することを示す。速度修正部136は、速度算出部130における演算の結果得られた粒子の速度を、エネルギ変化量算出部134によって算出された粒子のエネルギの変化量に基づいて変更する。特に速度修正部136は、速度算出部130における演算の結果得られた粒子の速度によって定まる運動エネルギにその粒子について算出されたエネルギの変化量を加えることで新たな運動エネルギを演算し、演算された新たな運動エネルギを与える粒子の速度を演算する。本実施の形態では、粒子系Sのi番目の粒子について、速度算出部130における演算の結果得られた粒子の速度viにより定まる運動エネルギにエネルギ変化量算出部134によって算出されたエネルギΔEiを加えることによってi番目の粒子の運動エネルギを修正する。そのように修正された運動エネルギに対応する速度を新たなi番目の粒子の速度とする。すなわち、速度修正部136は以下の式6によってi番目の粒子の速度を修正する。ここでviはi番目の粒子の修正された速度である。位置算出部138は、速度修正部136によって修正された粒子の速度に基づいて粒子の位置を算出する。位置算出部138は、粒子系Sのi番目の粒子について、所定の数値解析の手法に基づき時間刻みΔtを使用して離散化された粒子の位置と速度の関係式に、速度修正部136によって修正された粒子の速度viを適用することによって、粒子の位置xiベクトル量を演算する。この演算には以前のステップで演算された粒子の位置が使用される。位置更新部126は、記憶部160に記憶される粒子系Sの各粒子の位置を、位置算出部138によって算出された位置で更新する。また位置更新部126は、記憶部160に記憶される粒子系Sの各粒子の速度を、速度修正部136によって修正された速度で更新してもよい。終了条件判定部128は、繰り返し演算部120における繰り返し演算を終了すべきか否かを判定する。繰り返し演算を終了すべき終了条件は、例えば繰り返し演算が所定の回数行われたことや、粒子系Sが定常状態に達したことや、外部から終了の指示を受け付けたことである。終了条件判定部128は、終了条件が満たされる場合、繰り返し演算部120における繰り返し演算を終了させる。終了条件判定部128は、終了条件が満たされない場合、処理を力演算部122に戻す。すると力演算部122は、位置更新部126によって更新された粒子の位置で再び力を演算する。結果提示部150は、演算条件設定部110によって設定された粒子系Sの解析結果をユーザに提示する。結果提示部150は、物理量演算部152と、再スケーリング部154と、描画制御部156と、を含む。物理量演算部152は、繰り返し演算部120における繰り返し演算が終了した後、記憶部160に記憶される粒子系Sの情報に基づき粒子系Sの各種物理量、例えば温度や圧力や応力などを演算する。再スケーリング部154は、モデル生成部114において粒子系Sの粒子が繰り込まれた系の粒子とされた場合、物理量演算部152によって演算された物理量を繰り込まれる前の系の物理量に変換する。特に再スケーリング部154は、物理量演算部152によって演算された物理量に、物理量ごとに定まるスケーリング係数を乗じることによって繰り込まれる前の系の物理量を得る。物理量には応力や温度など繰り込み変換に際して不変となる物理量もあり、そのような物理量についてはスケーリング係数として1変換前後で不変が設定される。再スケーリング部154は、変換された物理量を出力装置104に表示させる。描画制御部156は、繰り返し演算部120における繰り返し演算が終了した後、記憶部160に記憶される粒子系Sの各粒子の位置、速度の情報に基づき、出力装置104に粒子系Sの時間発展や定常状態の様子をグラフィカルに表示させる。上述の実施の形態において、記憶部160の例は、ハードディスクやメモリである。また、本明細書の記載に基づき、各部を、図示しないCPUや、インストールされたアプリケーションプログラムのモジュールや、システムプログラムのモジュールや、ハードディスクから読み出したデータの内容を一時的に記憶するメモリなどにより実現できることは本明細書に触れた当業者には理解されるところである。以上の構成による解析装置100の動作を説明する。図2は、解析装置100における一連の処理の一例を示すフローチャートである。熱伝導率取得部112は、解析対象に対して設定されるべき熱伝導率を取得するS202。モデル生成部114は、解析対象を記述した粒子系Sを生成するS204。力演算部122は、粒子間の距離から粒子に働く力を演算するS206。速度算出部130は、演算された力を含む粒子の運動方程式から速度を演算するS208。熱流束密度算出部132は、取得された熱伝導率を使用して、粒子間の熱流束密度を演算するS210。エネルギ変化量算出部134は、演算された熱流束密度に起因する粒子のエネルギの変化量を演算するS212。速度修正部136は、運動方程式から演算された粒子の速度を、演算された粒子のエネルギの変化量を使用して修正するS214。位置算出部138は、修正された速度から粒子の位置を算出するS216。位置更新部126は、記憶部160に記憶される粒子の位置を算出された位置で更新するS218。終了条件判定部128は、終了条件が満たされるか否かを判定するS220。終了条件が満たされない場合S220のN、処理はS206に戻される。終了条件が満たされる場合S220のY、結果提示部150は演算結果をユーザに提示するS222。本実施の形態に係る解析装置100では、MD法やRMD法を使用して粒子系のふるまいを演算する際、隣接する粒子間に文献や実測等から得られる熱伝導率に基づく熱の流れを仮定し、演算に修正を加える。この熱の流れを演算する際の前提として、各粒子に温度が定義される。この手法では瞬間的に熱エネルギの流れが無秩序となる可能性があるが、解析装置100では時間刻みΔtごとに演算を繰り返すので、その繰り返しの回数を適宜増やすことで時間的なアンサンブルを十分に取り、信頼性の高い統計量を導くことができる。また本実施の形態に係る解析装置100によると、解析対象の熱伝導、特に従来のMD法やRMD法では取り扱いが困難な自由電子による熱伝導を補正しながらシミュレーションを行うことができる。その結果、相転移を伴わない通常の解析対象については熱伝導の側面からより正確にシミュレーションを行うことができる。また、一般にMD法やRMD法に基づく解析は相転移をより忠実に再現できるのであるが、本実施の形態に係る解析装置100によるとさらに熱伝導現象の再現性も向上するので、結果として液体金属の冷却凝固の過程や固体金属の液化の過程などの熱伝導がクリティカルとなる一次相転移の過程をより正確にシミュレーションできる。なお、粒子法において、粒子それぞれに温度をパラメータとして与え、熱伝導方程式を解く手法も考えられるが、この場合は潜熱などが取り扱えないので相転移を正確にシミュレートすることができない。これに対して本実施の形態に係る解析装置100では熱伝導現象および相転移の両方を好適に取り扱える。また本実施の形態に係る解析装置100によると、RMD法に倣って繰り込まれた粒子系が設定される場合、ミクロスケール、マクロスケールを問わず任意のスケールのシミュレーションが可能となる。また本実施の形態に係る解析装置100では、式3、式4、式5、式6など、物理的な根拠がある関係式にしたがって粒子の速度が修正される。これにより、特に熱に関して物理法則に則ったより忠実なシミュレーションが実現される。図3は、演算条件設定部110によって設定される粒子系10の一例を示す模式図である。x、y、zの領域は、それぞれ2.86×2.476×4.46nmである.図3のように、z方向の上端を100K、下端を200Kの温度で一定となるようにスケーリングを行う。このときスケーリング前後のエネルギ差から流入、流出エネルギを求める。図4は、粒子系10の下端からの流入エネルギ12および上端からの流出エネルギ14を時間の関数として示すグラフである。図4の縦軸はエネルギ、横軸は時間を示し、この場合の時間は繰り返し演算の回数に比例する。粒子系10が定常状態になると、図4に示されるように下端からの流入エネルギの絶対値と上端からの流出エネルギの絶対値とが同じになる。これらのエネルギが求まると、下記の式7から熱伝導率が求まる。図3、図4の例では、熱伝導率としてアルミニウムの238W/mKを使用した。図4から求めたエネルギの値を式7に代入し、熱伝導率を求めたところ、約238±20.0W/mKという結果を得たので、熱伝導率について自己無撞着selfconsistentであることが確認された。また、本実施の形態に係る手法を用いない、通常のMD法で熱伝導率を求めた場合、値は約3.3W/mKとなった。このことから、本実施の形態に係る手法は目的とする物質の熱伝導率を正しく再現できることが示された。以上、実施の形態に係る解析装置100の構成と動作について説明した。これらの実施の形態は例示であり、その各構成要素や各処理の組み合わせにいろいろな変形例が可能なこと、またそうした変形例も本発明の範囲にあることは当業者に理解されるところである。数値解析の手法にはVerlet法のように、粒子の位置を演算する際に粒子に働く力から粒子の位置を直接演算し、粒子の速度は陽に計算しなくてもよい手法がある。本実施の形態に係る技術的思想をこのような手法に適用する場合、例えば速度修正部は力に時間の次元を有する量を乗算した量に関して上記と同様の修正を行ってもよい。位置算出部は、そのように修正された量に基づいて粒子の位置を演算してもよい。