特許 6053418 解析方法および解析装置

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

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

(11)【特許番号】6053418

(24)【登録日】2016年12月9日

(45)【発行日】2016年12月27日

(54)【発明の名称】解析方法および解析装置

(51)【国際特許分類】

   G06F 19/00 20110101AFI20161219BHJP

【ＦＩ】

   G06F19/00 110

【請求項の数】6

【全頁数】14

(21)【出願番号】特願2012-208721(P2012-208721)

(22)【出願日】2012年9月21日

(65)【公開番号】特開2014-63388(P2014-63388A)

(43)【公開日】2014年4月10日

【審査請求日】2015年1月15日

(73)【特許権者】

【識別番号】000002107

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

(74)【代理人】

【識別番号】100105924

【弁理士】

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

(74)【代理人】

【識別番号】100109047

【弁理士】

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

(74)【代理人】

【識別番号】100109081

【弁理士】

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

(74)【代理人】

【識別番号】100116274

【弁理士】

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

(72)【発明者】

【氏名】大西 良孝

【審査官】 関 博文

(56)【参考文献】

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

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

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

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

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

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

【文献】 河野 瑛、田浦 健次朗，重心ボロノイ分割を用いた並列粒子法のための動的負荷分散法，情報処理学会研究報告 ２０１１（平成２３）年度▲２▼ ［ＣＤ−ＲＯＭ］，日本，一般社団法人情報処理学会，２０１１年 ８月１５日，Ｎｏ．１３０，（６６）

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

Ｇ０６Ｆ １９／００

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

【請求項1】

仮想空間に定義される粒子系を解析する解析方法であって、
当該解析方法は、解析装置において実行され、
前記解析装置が、粒子系に含まれる２つの粒子の間の面積を、２つの粒子の間の距離のみを引数とする関数に基づいて決定するステップと、
前記解析装置が、決定された面積を使用して粒子系の状態を更新するステップと、を含むことを特徴とする解析方法。

【請求項2】

前記決定するステップにおいて決定される面積は、２つの粒子の間の距離が大きいほど小さいことを特徴とする請求項１に記載の解析方法。

【請求項3】

前記決定するステップにおいて決定される面積は、２つの粒子の間の距離の非線形関数であることを特徴とする請求項１または２に記載の解析方法。

【請求項4】

前記更新するステップは、
ポテンシャルエネルギ関数を使用して粒子に働く力を演算する際、その粒子に力を及ぼす粒子をその粒子との距離が所定のカットオフ距離よりも小さな粒子に制限するステップと、
演算された粒子に働く力を含み粒子の運動を支配する支配方程式に基づき粒子の状態を更新するステップと、を含み、
前記決定するステップにおいて決定される面積は、２つの粒子の間の距離がカットオフ距離よりも大きい場合、実質的にゼロとなることを特徴とする請求項１から３のいずれかに記載の解析方法。

【請求項5】

仮想空間に定義される粒子系を解析する解析装置であって、
粒子系に含まれる２つの粒子の間の面積を、２つの粒子の間の距離のみを引数とする関数に基づいて決定する面積決定部と、
前記面積決定部によって決定された面積を使用して粒子系の状態を更新する状態更新部と、を備えることを特徴とする解析装置。

【請求項6】

仮想空間に定義される粒子系を解析する機能をコンピュータに実現させるコンピュータプログラムであって、
粒子系に含まれる２つの粒子の間の面積を、２つの粒子の間の距離のみを引数とする関数に基づいて決定する機能と、
決定された面積を使用して粒子系の状態を更新する機能と、を前記コンピュータに実現させることを特徴とするコンピュータプログラム。

【発明の詳細な説明】

【技術分野】

【0001】

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

【背景技術】

【0002】

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

【0003】

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

【0004】

そこで本出願人は特許文献１においてその解決法のひとつを提案している。

【先行技術文献】

【特許文献】

【0005】

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

【非特許文献】

【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】

以下、本実施の形態の原理を説明する。
粒子系の粒子を比較的密に詰めると、粒子系は面心立方格子（ｆｃｃ）を形成する傾向にある。したがって、ｆｃｃ構造を起点とする。
図１は、粒子系がｆｃｃ構造を有する場合のボロノイ面６を示す模式図である。ｆｃｃ構造を有する粒子系において、第１粒子２は第２粒子４の最近接粒子となっている。この粒子系を粒子の位置を母点としてボロノイ分割した場合、第１粒子２と第２粒子４との間にはボロノイ面６が定義される。このボロノイ面６の面積Ｓ_０は以下の式１により与えられる。

【数1】

ただし、ｒ_０は第１粒子２と第２粒子４との最安定距離である。なお、ボロノイ面６の面積Ｓ_０は粒子の断面積として扱われてもよい。

【0018】

第１粒子２を粒子系のｉ番目の粒子、第２粒子４を粒子系のｊ番目の粒子とするとき、本実施の形態では、粒子系の構造がｆｃｃ構造とは異なる場合でも、第１粒子２と第２粒子４との間の断面積ΔＳ_ｉｊを第１粒子２と第２粒子４との距離ｒ_ｉｊのみから求めることができると仮定する。すなわち、断面積ΔＳ_ｉｊは第１粒子２、第２粒子４以外の他の粒子の位置に依存しないと仮定する。

【0019】

図２は、第１粒子２と第２粒子４との間の断面積ΔＳ_ｉｊとそれらの粒子の間の距離ｒ_ｉｊとの間に仮定される関係の一例を示すグラフである。大局的には、断面積ΔＳ_ｉｊは粒子間距離ｒ_ｉｊが大きいほど小さくなるよう設定される。断面積ΔＳ_ｉｊは粒子間距離ｒ_ｉｊの線形関数８とされてもよい。または、断面積ΔＳ_ｉｊは粒子間距離ｒ_ｉｊの非線形関数とされてもよく、特に２次以上の関数１０とされてもよい。また、断面積ΔＳ_ｉｊは粒子間距離ｒ_ｉｊが後述のカットオフ距離ｒ_ｃよりも大きい場合、実質的にゼロとなるよう設定される。

【0020】

本実施の形態に係る解析装置では、ユーザは、粒子間の断面積と粒子間距離とに図２に示されるような１対１の関係を仮定し、仮定された関係すなわち関数を解析装置に登録する。解析装置は、数値演算の過程で粒子系に含まれる２つの粒子の間の断面積を求める必要がある場合、登録された関係を使用することで、その断面積を他の粒子の位置に依らずに決定する。解析装置は、決定された断面積を使用して数値演算を継続し、粒子系の状態を更新する。これにより、２つの粒子の間の断面積を求める必要がある場合に毎回比較的厳密にその断面積を導出する場合と比較して、粒子間の断面積の導出に係る計算負荷を低減することができる。その結果、ユーザはより速く解析結果を得ることができる。

【0021】

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

【0022】

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

【0023】

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

【0024】

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

【0025】

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

【0026】

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

【0027】

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

【0028】

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

【0029】

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

【0030】

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

【0031】

面積決定部１１２は、粒子系に含まれる２つの粒子の間の断面積を、他の粒子の位置に依らずに決定する。特に面積決定部１１２は、粒子データ保持部１１４によって保持される粒子系のデータを参照し、粒子系のｉ番目の粒子とｊ番目の粒子と（１≦ｉ、ｊ≦Ｎ）の粒子間距離ｒ_ｉｊを演算する。面積決定部１１２は、図２に示される２次以上の関数１０に演算された粒子間距離ｒ_ｉｊを代入することにより、ｉ番目の粒子とｊ番目の粒子との間の断面積ΔＳ_ｉｊを演算する。

【0032】

熱伝導演算部１１６は、離散化された熱伝導方程式に基づいて各粒子の温度を演算する。熱伝導演算部１１６では、有限体積法を用いた温度場の解析が実行されてもよい。
以下、熱伝導演算部１１６で使用される離散化された熱伝導方程式の導出過程を示す。
熱伝導方程式は以下の微分方程式（式２）で与えられる。

【数2】

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

【0033】

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

【数3】

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

【数4】

となる。

【0034】

式４を離散化すると、

【数5】

となる。ここで、ΔＶ_ｉはｉ番目の粒子の位置を中心とする半径ｒ_０の球の体積（４πｒ_０^３／３）、ΔＳ_ｉｊは面積決定部１１２によって決定されるｉ番目の粒子とｊ番目の粒子との間の断面積、ｒ_ｉｊはｉ番目の粒子とｊ番目の粒子との距離である。また、Ｔ_ｉ^ｎは繰り返し演算のｎ回目におけるｉ番目の粒子の温度、Ｋ_ｉｊはｉ番目の粒子とｊ番目の粒子との間の熱伝導率である。熱伝導率Ｋ_ｉｊは、面積決定部１１２における断面積の近似処理に伴い文献値から補正される。この補正で使用される係数は粒子系の構造の状態に依存し、０．２〜２．０ファクター倍程度とすると大抵の場合理論値と良い一致を示す。

【0035】

式５は、熱伝導演算部１１６で使用される離散化された熱伝導方程式である。式５により、面積決定部１１２によって決定される断面積と、ｎ回目の演算により決定された各粒子の温度と、粒子間距離と、から、ｎ＋１回目の演算におけるｉ番目の粒子の温度を決定することができる。なお、Ｔ_ｉ^０には、温度付与部１３４によって付与された初期温度を使用する。

【0036】

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

【0037】

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

【0038】

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

【0039】

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

【0040】

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

【数6】

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

【表1】

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

【0041】

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

【数7】

の標準偏差σを有する。

【0042】

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

【数8】

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

【0043】

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

【0044】

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

【0045】

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

【0046】

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

【0047】

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

【0048】

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

【0049】

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

【0050】

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

【0051】

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

【0052】

本実施の形態に係る解析装置１００によると、離散化された熱伝導方程式を演算する際の２粒子間の断面積ΔＳ_ｉｊは、その２粒子間の距離ｒ_ｉｊに基づき他の粒子の位置に依らずに決定される。したがって、例えばボロノイ解析等のより厳密な手法により断面積を求める場合と比較して、計算の負荷をかなり低減することができる。

【0053】

一般に粒子系の構造がｆｃｃ構造から離れるほどボロノイ解析による計算負荷は大きくなるが、本実施の形態に係る手法によって断面積を決定する際の計算負荷は基本的には粒子系の構造には依存しない。したがって、本実施の形態に係る手法は、粒子系の構造がｆｃｃ構造からかなり離れる場合や、粒子系の構造が時間と共に大きく変化する場合により好適に適用されうる。より具体的には、解析対象が流体的な振る舞いを示す場合に、本実施の形態に係る解析手法が好適に適用されうる。

【0054】

定量的で正確な解析結果を得たい場合などには、ボロノイ解析等を使用して断面積を決定するほうが好ましい。しかしながら、解析装置１００による解析結果の利用方法としては、例えば定性的な議論をする際の参考にする場合など、それほど正確な結果は必要とされずむしろより早く結果を得たい場合も多くある。そこで、本実施の形態に係る解析装置１００は、ある程度厳密性を犠牲にしてその分より早く解析結果を提供することで、そのような要望に対応できる。

【0055】

また、本実施の形態に係る手法では、断面積は粒子間距離が大きくなるほど小さくなるよう設定され、特に粒子間距離がカットオフ距離ｒ_ｃを超えると断面積は実質的にゼロになる。これは、比較的遠い粒子の間での熱のやり取りは小さい、という物理的見地に合致する仮定である。また、ある粒子と相互作用する粒子をカットオフ距離ｒ_ｃ内の粒子に制限しているので、ある粒子と熱をやり取りする粒子を同じくカットオフ距離ｒ_ｃ内の粒子に制限することもまた、物理的見地に合致する仮定である。このように物理的見地に合致する仮定の下で得られる解析結果は、より高い物理的な整合性を有することが期待される。

【0056】

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

【0057】

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

【0058】

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

【0059】

図６は、確認計算としての非定常解析で使用された粒子系３００を示す模式図である。粒子系３００はアモルファス金属の棒を模している。アモルファス金属では一般に自由電子からの熱伝導への寄与が無視できず、また結晶構造も流体的に変化する。この棒の両端の温度は０（Ｋ）で固定され、初期の温度分布は以下の式９により与えられる。

【数9】

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

【0060】

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

【数10】

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

【数11】

【0061】

図７は、本実施の形態に係る手法を使用した計算結果を示すグラフである。本実施の形態に係る手法によると、粒子系３００が時間発展を始めてから１（μｓ）経過後、２（μｓ）経過後、３（μｓ）経過後および４（μｓ）経過後のいずれにおいても、温度分布の計算値（点で表示）と理論値（実線で表示）とが良く一致していることが分かる。

【0062】

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

【0063】

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

【符号の説明】

【0064】

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

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】