特許 6129193 解析装置

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】6129193

(24)【登録日】2017年4月21日

(45)【発行日】2017年5月17日

(54)【発明の名称】解析装置

(51)【国際特許分類】

   G06F 19/00 20110101AFI20170508BHJP

【ＦＩ】

   G06F19/00 110

【請求項の数】4

【全頁数】14

(21)【出願番号】特願2014-536606(P2014-536606)

(86)(22)【出願日】2013年9月20日

(86)【国際出願番号】JP2013005581

(87)【国際公開番号】WO2014045594

(87)【国際公開日】20140327

【審査請求日】2015年11月17日

(31)【優先権主張番号】特願2012-208720(P2012-208720)

(32)【優先日】2012年9月21日

(33)【優先権主張国】JP

(73)【特許権者】

【識別番号】000002107

【氏名又は名称】住友重機械工業株式会社

(74)【代理人】

【識別番号】100105924

【弁理士】

【氏名又は名称】森下 賢樹

(74)【代理人】

【識別番号】100109047

【弁理士】

【氏名又は名称】村田 雄祐

(74)【代理人】

【識別番号】100109081

【弁理士】

【氏名又は名称】三木 友由

(74)【代理人】

【識別番号】100116274

【弁理士】

【氏名又は名称】富所 輝観夫

(72)【発明者】

【氏名】大西 良孝

(72)【発明者】

【氏名】市嶋 大路

【審査官】 松野 広一

(56)【参考文献】

【文献】 特開２０１１−１７５３２７（ＪＰ，Ａ）

【文献】 特開２０１０−１７０３０９（ＪＰ，Ａ）

【文献】 特開２０１２−１２８４９０（ＪＰ，Ａ）

【文献】 特開２０１２−１５０６７３（ＪＰ，Ａ）

【文献】 特開２０１０−２３０３３１（ＪＰ，Ａ）

【文献】 特開２００１−３３４３０５（ＪＰ，Ａ）

【文献】 特開昭６１−１３２２０５（ＪＰ，Ａ）

(58)【調査した分野】（Int.Cl.，ＤＢ名）

Ｇ０６Ｆ １９／００

Ｇ０１Ｎ １５／００

(57)【特許請求の範囲】

【請求項1】

仮想空間に定義される粒子系の各粒子の運動を支配する支配方程式を数値的に演算することによりその粒子系を解析する解析装置であって、
粒子系の各粒子のパラメータのひとつである温度を演算する温度演算部と、
各粒子をその粒子について前記温度演算部によって演算された温度の熱浴に浸すと仮定した場合の、その粒子に働く力を演算する力演算部と、
前記力演算部によって演算された力に基づき各粒子の状態を更新する状態更新部と、を備え、
前記温度演算部は、
各粒子の位置に基づいて仮想空間にボロノイ多面体を生成するボロノイ処理部と、
前記ボロノイ処理部によって生成されたボロノイ多面体を単位として離散化された熱伝導方程式に基づいて各粒子の温度を演算する熱伝導演算部と、を含み、
力演算部は、
粒子間の距離に基づき各粒子に働く力を演算する粒子間作用演算部と、
前記粒子間作用演算部によって演算された力に粘性力およびランダム力を付加する熱浴作用演算部と、を含むことを特徴とする解析装置。

【請求項2】

前記熱伝導方程式は、粒子系の変形による発熱を表す項を含むことを特徴とする請求項１に記載の解析装置。

【請求項3】

粒子系の粒子が金属の粒子を模するように各粒子についての物質定数が設定されることを特徴とする請求項１または２に記載の解析装置。

【請求項4】

仮想空間に定義される粒子系の各粒子の運動を支配する支配方程式を数値的に演算することによりその粒子系を解析する機能をコンピュータに実現させるコンピュータプログラムであって、
粒子系の各粒子のパラメータのひとつである温度を演算する温度演算機能と、
各粒子をその粒子について演算された温度の熱浴に浸すと仮定した場合の、その粒子に働く力を演算する力演算機能と、
演算された力に基づき各粒子の状態を更新する状態更新機能と、を前記コンピュータに実現させ、
前記温度演算機能は、
各粒子の位置に基づいて仮想空間にボロノイ多面体を生成するボロノイ処理機能と、
前記ボロノイ処理機能によって生成されたボロノイ多面体を単位として離散化された熱伝導方程式に基づいて各粒子の温度を演算する熱伝導演算機能と、を含み、
前記力演算機能は、
粒子間の距離に基づき各粒子に働く力を演算する粒子間作用演算機能と、
前記粒子間作用演算機能によって演算された力に粘性力およびランダム力を付加する熱浴作用演算機能と、を含むことを特徴とするコンピュータプログラム。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、粒子系を解析する解析装置に関する。

【背景技術】

【0002】

古典力学や量子力学等を基に計算機を用いて物質科学全般の現象を探るための方法として、分子動力学法（Molecular Dynamics Method、以下ＭＤ法と称す）や、量子分子動力学法（Quantum Molecular Dynamics Method）や、ＭＤ法をマクロスケールの系を扱えるように発展させた繰り込み群分子動力学法（Renormalized Molecular Dynamics、以下ＲＭＤ法と称す）に基づくシミュレーションが知られている。

【0003】

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

【0004】

そこで本出願人は特許文献１および特許文献２において、いくつかの解決法を提案している。

【先行技術文献】

【特許文献】

【0005】

【特許文献1】特開２０１０−１７０３０９号公報

【特許文献2】特開２０１２−１５０６７３号公報

【非特許文献】

【0006】

【非特許文献1】John C. Tully、「Dynamics of gas-surface interaction: 3D generalized Langevin model applied to fcc and bcc surface」、J. Chem. Phys.、1975、73

【発明の概要】

【発明が解決しようとする課題】

【0007】

特許文献１に記載の手法では、粒子に温度のパラメータを持たせ、有限体積法（Finite Volume Method、ＦＶＭ）により熱伝導方程式を解くことで温度場を求めている。しかしながら、その温度場は粒子の速度の分散を反映していない場合が多く、適用範囲は比較的限られている。

【0008】

特許文献２には、粒子の運動エネルギを確率変数に換算し、その勾配から、フーリエ則に従ってエネルギを受け渡す技術が開示されている。しかしながら、温度の勾配によってはエネルギが大きくなり、粒子の持つエネルギがマイナスになる場合がある。それを防ぐには時間刻みを小さく取ればよいのだが、そうすると計算負荷が増大しうる。

【0009】

本発明はこうした課題に鑑みてなされたものであり、その目的は、複数の粒子を含む粒子系を使用したシミュレーションにおいて、熱伝導現象をより正確に扱えるようにする解析技術の提供にある。

【課題を解決するための手段】

【0010】

本発明のある態様は、解析装置に関する。この解析装置は、仮想空間に定義される粒子系の粒子の運動を支配する支配方程式を数値的に演算することによりその粒子系を解析する解析装置であって、粒子系の粒子のパラメータのひとつである温度を演算する温度演算部と、粒子をその粒子について温度演算部によって演算された温度の熱浴に浸すと仮定した場合の、その粒子に働く力を演算する力演算部と、力演算部によって演算された力に基づき粒子の状態を更新する状態更新部と、を備える。

【0011】

この態様によると、粒子系の熱伝導をより好適に扱うことができる。

【0012】

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

【発明の効果】

【0013】

本発明によれば、複数の粒子を含む粒子系を使用したシミュレーションにおいて、熱伝導現象をより正確に扱うことができる。

【図面の簡単な説明】

【0014】

【図1】第１の実施の形態に係る解析装置の機能および構成を示すブロック図である。

【図2】図１の粒子データ保持部の一例を示すデータ構造図である。

【図3】図１の解析装置における一連の処理の一例を示すフローチャートである。

【図4】第１の実施の形態に係る手法の確認計算で使用された粒子系を示す模式図である。

【図5】第１の実施の形態に係る手法を使用した計算結果を示すグラフである。

【図6】第１の実施の形態に係る手法を使用しない場合の計算結果を示すグラフである。

【図7】図７（ａ）〜（ｆ）は、第２の実施の形態に係る手法を使用した計算結果の時間変化を示す模式図である。

【発明を実施するための形態】

【0015】

以下、各図面に示される同一または同等の構成要素、部材、処理には、同一の符号を付するものとし、適宜重複した説明は省略する。

【0016】

実施の形態に係る解析装置は、解析対象を複数の粒子を含む粒子系で記述し、粒子の運動方程式を数値的に演算することにより粒子系を解析する。その際、解析装置は熱伝導方程式を解いて温度場を求める。解析装置は、粒子の運動方程式に、求められた温度と実質的に等しい温度の熱浴に粒子が浸されると仮定したときに得られる粒子に働く力を導入する。具体的には、熱浴の粒子からのランダム力と、粘性力によるダンパーと、を与えるＬａｎｇｅｖｉｎ法（非特許文献１参照）が適用される。これにより、自由電子からの熱伝導への寄与が無視できない金属などの解析対象の熱伝導を再現する際の物理的な整合性を高めることができる。

【0017】

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

【0018】

本実施の形態ではＭＤ法に倣って粒子系を解析する場合について説明するが、ＲＭＤ法に倣って粒子系を解析する場合や、ＤＥＭ（Distinct Element Method）やＳＰＨ（Smoothed Particle Hydrodynamics）やＭＰＳ（Moving Particle Semi-implicit）などの粒子法に倣って粒子系を解析する場合にも、本実施の形態に係る技術的思想を適用できることは本明細書に触れた当業者には明らかである。

【0019】

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

【0020】

解析装置１００は、粒子系取得部１０８と、温度付与部１３４と、繰り返し演算部１２０と、表示制御部１１８と、粒子データ保持部１１４と、を備える。

【0021】

粒子系取得部１０８は、入力装置１０２を介してユーザから取得する入力情報に基づき、１、２または３次元の仮想空間内に定義されるＮ（Ｎは自然数）個の粒子からなる粒子系のデータを取得する。粒子系の粒子は分子または原子に対応してもよい。

【0022】

粒子系取得部１０８は、入力情報に基づき仮想空間内にＮ個の粒子を配置し、配置されたそれぞれの粒子に速度を付与する。すなわち、粒子系取得部１０８は、粒子系に初期位置および初期速度を付与する。粒子系取得部１０８は、配置された粒子を特定する粒子ＩＤと、その粒子の位置と、その粒子の速度と、を対応付けて粒子データ保持部１１４に登録する。

【0023】

温度付与部１３４は、入力装置１０２を介してユーザから取得する入力情報に基づき、粒子系取得部１０８によって取得された粒子系の粒子に温度を付与する。このように付与される温度は粒子のパラメータのひとつである。例えば、温度付与部１３４は、ディスプレイ１０４を介してユーザに、粒子系の粒子の温度の初期値の入力を求める。温度付与部１３４は、入力された温度の初期値を粒子ＩＤと対応付けて粒子データ保持部１１４に登録する。

【0024】

以下では粒子系の粒子は全て同質または同等なものとして設定され、かつ、ポテンシャルエネルギ関数は２体のポテンシャルであって粒子によらずに同じ形を有するものとして設定される場合について説明する。しかしながら、他の場合にも本実施の形態に係る技術的思想を適用できることは、本明細書に触れた当業者には明らかである。

【0025】

繰り返し演算部１２０は、粒子データ保持部１１４によって保持されるデータが表す粒子系の各粒子の運動を支配する支配方程式を数値的に演算する。特に繰り返し演算部１２０は、離散化された粒子の運動方程式にしたがった繰り返し演算を行う。
繰り返し演算部１２０は、温度演算部１１０と、力演算部１２２と、粒子状態演算部１２４と、状態更新部１２６と、終了条件判定部１２８と、を含む。

【0026】

温度演算部１１０は、連続体近似を用いて粒子系の各粒子の温度を演算する。特に温度演算部１１０は、離散化された熱伝導方程式に基づき粒子の温度を演算する。温度演算部１１０は、ボロノイ分割部１１２と、熱伝導演算部１１６と、を有する。

【0027】

ボロノイ分割部１１２は、粒子系が定義されている仮想空間の部分をボロノイ分割する。すなわち、ボロノイ分割部１１２は、各粒子の位置に基づいて仮想空間内にボロノイ多面体を生成する。ボロノイ分割部１１２はまず、各粒子を頂点として仮想空間の部分を３次元デローニ（ドロネー）分割する。次にボロノイ分割部１１２は、デローニ分割の結果得られる四面体要素からボロノイ多面体要素を生成する。これにより、仮想空間の部分は粒子の位置を母点としてボロノイ分割される。

【0028】

熱伝導演算部１１６は、ボロノイ分割部１１２によって生成されたボロノイ多面体を単位として離散化された熱伝導方程式に基づいて各粒子の温度を演算する。熱伝導演算部１１６では、有限体積法を用いた温度場の解析が実行されてもよい。熱伝導演算部１１６は、温度場の解析において、ボロノイ分割部１１２によって生成されたボロノイ多面体の面積や体積、および所定の微小な時間刻みΔｔだけ前の時刻の温度場の情報（すなわち、繰り返し演算における１サイクル前の温度場の情報）を使用する。特に熱伝導演算部１１６は、ボロノイ多面体のそれぞれをＦＶＭのひとつのコントロールボリュームとして解析する。

【0029】

以下、熱伝導演算部１１６で使用される離散化された熱伝導方程式の導出過程を示す。
熱伝導方程式は以下の微分方程式（式１）で与えられる。

【数1】

ここで、ρは密度、Ｃｖは比熱、Ｔは温度、ｔは時間、λは熱伝導率、Ｑは単位体積当たりの発熱量である。

【0030】

式１の両辺を体積積分すると、

【数2】

となる。式２の右辺第１項にガウスの定理を適用すると、

【数3】

となる。

【0031】

式３をボロノイ多面体を単位として離散化すると、

【数4】

となる。ここで、ΔＶ_ｉはｉ番目の粒子の位置を母点とするボロノイ多面体の体積、ΔＳ_ｉｊはｉ番目の粒子とｊ番目の粒子との間にあるボロノイ面の面積、ｒ_ｉｊはｉ番目の粒子とｊ番目の粒子との距離である。また、Ｔ_ｉ^ｎは時間に関する繰り返し演算のｎ回目におけるｉ番目の粒子の温度（ｉ番目の粒子が属するボロノイ多面体の温度とも言える）、λ_ｉｊはｉ番目の粒子とｊ番目の粒子との間の熱伝導率である。式４は、熱伝導演算部１１６で使用される離散化された熱伝導方程式である。式４により、ボロノイ多面体の体積および面積と、ｎ回目の演算により決定された各粒子の温度と、粒子間距離と、から、ｎ＋１回目の演算におけるｉ番目の粒子の温度を決定することができる。なお、Ｔ_ｉ^０には、温度付与部１３４によって付与された初期温度を使用する。

【0032】

力演算部１２２は、粒子をその粒子について温度演算部１１０によって演算された温度の熱浴に浸すと仮定した場合の、その粒子に働く力を演算する。力演算部１２２は、粒子間作用演算部１３０と、熱浴作用演算部１３２と、を有する。

【0033】

粒子間作用演算部１３０は、粒子データ保持部１１４によって保持される粒子系のデータを参照し、粒子系の各粒子について、粒子間の距離に基づきその粒子に働く力を演算する。粒子間作用演算部１３０は、粒子系のｉ番目（１≦ｉ≦Ｎ）の粒子について、そのｉ番目の粒子との距離が所定のカットオフ距離よりも小さな粒子（以下、近接粒子と称す）を決定する。

【0034】

粒子間作用演算部１３０は、各近接粒子について、その近接粒子とｉ番目の粒子との間のポテンシャルエネルギ関数およびその近接粒子とｉ番目の粒子との距離に基づいて、その近接粒子がｉ番目の粒子に及ぼす力を演算する。特に粒子間作用演算部１３０は、その近接粒子とｉ番目の粒子との距離の値におけるポテンシャルエネルギ関数のグラジエント（Gradient）の値から力を算出する。粒子間作用演算部１３０は、近接粒子がｉ番目の粒子に及ぼす力を全ての近接粒子について足し合わせることによって、ｉ番目の粒子に働く力を算出する。粒子間作用演算部１３０によって算出される力は、粒子間の相互作用に基づく力である。

【0035】

ｉ番目の粒子について温度演算部１１０によって演算された温度Ｔ_ｉで実質的に定温の熱浴にｉ番目の粒子を浸すと仮定した場合、Ｌａｎｇｅｖｉｎ法によるとｉ番目の粒子に以下の２つの力が働く。

【0036】

（１）ダンパーによる力（粘性力）
熱浴の粘性が粒子に作用する力である。粘性力の減衰定数αは、デバイ周波数ω_Ｄと粒子の質量ｍとを用いて、

【数5】

と表される。以下の表は、代表的な金属物質の減衰定数αの値を示す。この表では粒子を原子に対応付けている。

【表1】

この表に示される通り、粒子が金属原子に対応する場合、減衰定数αは１．０ｘ１０^−１２（ｋｇ／ｓ）のオーダーとなる。なお、デバイ周波数ω_Ｄは粒子の質量に依存するので、粒子の質量が原子の質量のβ倍になると、減衰定数αはβ^０．５倍となる。例えば、鉄の物性を持ちながら質量が鉄原子の１００倍大きい粒子を用いる場合、表１に従うと減衰定数αは２．９９ｘ１０^−１１（ｋｇ／ｓ）となる。

【0037】

（２）ランダム力
熱浴中の粒子が衝突する力に相当する。ランダム力は、

【数6】

の標準偏差σを有する。ここで、Ｋｂはボルツマン定数である。

【0038】

熱浴作用演算部１３２は、粒子系のｉ番目（１≦ｉ≦Ｎ）の粒子について、式５および式６に基づき粘性力とランダム力とを演算する。熱浴作用演算部１３２は、粒子間の相互作用に基づくｉ番目の粒子に働く力に、ｉ番目の粒子について演算された粘性力およびランダム力を加えることによって、ｉ番目の粒子に働くトータルの力を演算する。ｉ番目の粒子に働くトータルの力Ｆ_ｉは以下の式７で表される。

【数7】

ここで、φ_ｉｊはｉ番目の粒子とｊ番目の粒子との間のポテンシャルエネルギ関数、ｖ_ｉはｉ番目の粒子の速度、Ｆ_{ｒａｎｄｏｍ}は標準偏差σを有するランダム力である。記号の上の矢印はベクトル量であることを示す。

【0039】

粒子状態演算部１２４は粒子データ保持部１１４に保持される粒子系のデータを参照し、粒子系の各粒子について、離散化された粒子の運動方程式に熱浴作用演算部１３２によって演算されたトータルの力を適用することによって粒子の位置および速度のうちの少なくともひとつを演算する。本実施の形態では、粒子状態演算部１２４は粒子の位置および速度の両方を演算する。

【0040】

粒子状態演算部１２４は、熱浴作用演算部１３２によって演算されたトータルの力を含む離散化された粒子の運動方程式から粒子の速度を演算する。粒子状態演算部１２４は、粒子系のｉ番目の粒子について、蛙跳び法やオイラー法などの所定の数値解析の手法に基づき時間刻みΔｔを使用して離散化された粒子の運動方程式に、ｉ番目の粒子について熱浴作用演算部１３２によって演算されたトータルの力を代入することによって、粒子の速度を演算する。この演算には以前の繰り返し演算のサイクルで演算された粒子の速度が使用される。

【0041】

粒子状態演算部１２４は、演算された粒子の速度に基づいて粒子の位置を算出する。粒子状態演算部１２４は、粒子系のｉ番目の粒子について、所定の数値解析の手法に基づき時間刻みΔｔを使用して離散化された粒子の位置と速度の関係式に、演算された粒子の速度を適用することによって、粒子の位置を演算する。この演算には以前の繰り返し演算のサイクルで演算された粒子の位置が使用される。

【0042】

状態更新部１２６は、粒子状態演算部１２４における演算結果に基づき粒子系の各粒子の状態を更新する。状態更新部１２６は、粒子データ保持部１１４に保持される粒子系の各粒子の位置および速度のそれぞれを、粒子状態演算部１２４によって演算された位置および速度で更新する。

【0043】

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

【0044】

表示制御部１１８は、粒子データ保持部１１４に保持されるデータが表す粒子系の各粒子の位置、速度、温度に基づき、ディスプレイ１０４に粒子系の時間発展の様子やある時刻における状態を表示させる。この表示は、静止画または動画の形式で行われてもよい。

【0045】

図２は、粒子データ保持部１１４の一例を示すデータ構造図である。粒子データ保持部１１４は、粒子ＩＤと、粒子の位置と、粒子の速度と、粒子の温度と、を対応付けて保持する。

【0046】

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

【0047】

以上の構成による解析装置１００の動作を説明する。
図３は、解析装置１００における一連の処理の一例を示すフローチャートである。解析装置１００は、粒子系の初期状態すなわち各粒子の初期位置、初期速度および初期温度を決定する（Ｓ２０２）。解析装置１００は、各粒子の位置に基づきボロノイ解析を実行し、ボロノイ多面体を生成する（Ｓ２０４）。解析装置１００は、ＦＶＭを使用して温度場を解析し（Ｓ２０６）、各粒子の温度を更新する。解析装置１００は、粒子間のポテンシャルエネルギ関数に基づく各粒子に働く力を演算する（Ｓ２０８）。解析装置１００は、ステップＳ２０８で演算された力に、粘性力およびランダム力を付加する（Ｓ２１０）。解析装置１００は、ステップＳ２１０で演算された力を含む粒子の運動方程式から粒子の速度と位置を演算する（Ｓ２１２）。解析装置１００は、粒子データ保持部１１４に保持される粒子の位置、速度を演算された位置、速度で更新する（Ｓ２１４）。解析装置１００は、終了条件が満たされるか否かを判定する（Ｓ２１６）。終了条件が満たされない場合（Ｓ２１６のＮ）、処理はステップＳ２０４に戻される。終了条件が満たされる場合（Ｓ２１６のＹ）、処理は終了する。

【0048】

温度演算部１１０では連続体近似により各粒子の温度が導出されるので、温度演算部１１０によって演算される粒子の温度は、本来の温度の定義である粒子の速度の分散と大きく異なることが多い。本発明者はそのような非整合性を軽減または除去するために、温度演算部１１０によって温度を求めた後、その温度に由来する運動エネルギを粒子の運動に反映させることに想到した。ここで、例えば温度スケーリングなどで粒子の速度を強制的に温度に対応する速度に変更することが考えられるが、これは運動を強制的に拘束するため、非物理的である。

【0049】

そこで、本実施の形態に係る解析装置１００では、粒子を温度演算部１１０によって演算された温度の熱浴に浸すと仮定することで、運動方程式の力の項を温度に基づき修正している。これにより、温度演算部１１０によって演算された温度を粒子の速度場に反映させることができ、温度演算部１１０によって演算された温度場により自然な形で導くことができる。その結果、物理的な矛盾がより少ないモデルを提供できる。

【0050】

本発明者は、本実施の形態に係る手法の確認計算を行った。本実施の形態が使用するＭＤ法では、基本的に熱伝導として粒子の格子振動のみしか扱えないため、自由電子の寄与が反映されない。したがって、ＭＤ法で金属の解析対象を扱う場合、すなわち粒子系の粒子が金属の粒子を模するように粒子についての物質定数（例えば、デバイ温度やデバイ周波数や原子量、および、熱伝導方程式における密度や比熱や熱伝導率）が設定される場合、本実施の形態に係る手法が非常に有効である。

【0051】

図４は、確認計算で使用された粒子系３００を示す模式図である。粒子系３００は金属の棒を模している。この棒の両端の温度は０（Ｋ）で固定され、初期の温度分布は以下の式８により与えられる。

【数8】

ここで、Ｌは棒の長さ、ｒは端からの距離、Ｔ（ｒ）は距離ｒにおける温度である。

【0052】

この場合、時間ｔが経過すると、温度分布の理論式は以下の式９により与えられる。

【数9】

ここで、ａは熱拡散定数であり、以下の式１０で与えられる関係が成り立つ。

【数10】

【0053】

図５は、本実施の形態に係る手法を使用した計算結果を示すグラフである。図６は、本実施の形態に係る手法を使用しない場合の計算結果を示すグラフである。本実施の形態に係る手法によると、粒子系３００が時間発展を始めてから０．３（ｎｓ）経過後、０．６（ｎｓ）経過後、および０．９（ｎｓ）経過後のいずれにおいても、温度分布の計算値（黒点で表示）と理論値（実線で表示）とが良く一致していることが分かる。これに対して本実施の形態に係る手法を使用しない通常のＭＤ法の場合、理論値と比較して熱の拡散が非常に遅くなっていることが分かる。

【0054】

（第２の実施の形態）
第２の実施の形態では、粒子系の構造の変化、すなわち粒子配置の変化に伴う発熱を考慮した場合について説明する。第２の実施の形態に係る解析装置１００は、図１と同様の構成を有する。以下、第１の実施の形態との相違点を中心に説明する。

【0055】

熱伝導演算部１１６は、第１の実施の形態と同様に、ボロノイ分割部１１２によって生成されたボロノイ多面体を単位として離散化された熱伝導方程式に基づいて各粒子の温度を演算する。熱伝導演算部１１６で使用される離散化された熱伝導方程式は、第１の実施の形態と同様に、式４で表される。本実施の形態では、粒子系の構造の変化に伴う発熱を考慮するため、式４の発熱量Ｑは以下の式１１で表される。

【数11】

ここで、φは粒子間のポテンシャルエネルギ関数、Ｋ_ｉはｉ番目の粒子の運動エネルギ、ｎは時間に関する繰り返し演算の回数、Δｔは時間刻み、Ｆ_ｒｉｊはｉ番目の粒子とｊ番目の粒子との間で生じる摩擦力、ｖ_ｒｉｊはｉ番目の粒子とｊ番目の粒子との相対速度である。式１１は、粒子系の変形により粒子のもつ全エネルギが増加すれば発熱することを示している。式４および式１１により、ｎ＋１回目の演算におけるｉ番目の粒子の温度を決定することができる。

【0056】

本実施の形態に係る解析装置１００によれば、第１の実施の形態と同様に、温度演算部１１０によって演算された温度を粒子の速度場に反映させることができ、温度演算部１１０によって演算された温度場により自然な形で導くことができる。その結果、物理的な矛盾がより少ないモデルを提供できる。加えて、本実施の形態に係る解析装置１００によれば、温度演算部１１０（熱伝導演算部１１６）により演算される温度に、粒子系の構造の変化による発熱を反映させることができる。そのため、粒子系の構造に変形が生じる場合についても、物理的な矛盾がより少ないモデルを提供できる。これにより、例えば塑性加工など金属変形に伴う発熱が関与する現象をより正確にシミュレートすることができ、加工時の温度上昇を予測することができる。

【0057】

本発明者は、本実施の形態に係る手法の確認計算を行った。図７（ａ）〜（ｆ）は、本実施の形態に係る手法を使用した計算結果を示す模式図である。粒子系４００は、０．６ｍｍ（Ｘ方向）×０．６ｍｍ（Ｙ方向）×０．９５ｍｍ（Ｚ方向）の金属の角材を模している。図７（ａ）〜（ｆ）は、この粒子系４００の上下１層にそれぞれ応力５０ＧＰａ相当の引っ張り力を加えた場合の時間変化を示している。図７（ａ）は、粒子系４００に引っ張り力を加えた瞬間を示している。図７（ｂ）、（ｃ）、（ｄ）、（ｅ）、（ｆ）はそれぞれ、粒子系４００が時間発展を始めてから２５μｓ経過後、５０μｓ経過後、７５μｓ経過後、１００μｓ経過後、１２５μｓ経過後の状態を示す。この結果から、引っ張り力が加えられ、粒子系４００の構造が変化したことにより、粒子系４００の温度が上昇していることが分かる。これは、粒子系４００の構造が変化すると熱が発生するという知見と合致する。

【0058】

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

【0059】

実施の形態では、繰り返し演算部１２０において粒子の位置と速度の両方を演算する場合について説明したが、これに限られない。例えば、数値解析の手法にはＶｅｒｌｅｔ法のように、粒子の位置を演算する際に粒子に働く力から粒子の位置を直接演算し、粒子の速度は陽に計算しなくてもよい手法もあり、本実施の形態に係る技術的思想をそのような手法に適用してもよい。

【符号の説明】

【0060】

１００ 解析装置、 １０２ 入力装置、 １０４ ディスプレイ、 １０８ 粒子系取得部、 １１４ 粒子データ保持部、 １１８ 表示制御部、 １２０ 繰り返し演算部。

【産業上の利用可能性】

【0061】

本発明によれば、複数の粒子を含む粒子系を使用したシミュレーションにおいて、熱伝導現象をより正確に扱うことができる。

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】