JP6472346B2 - 粉粒体のシミュレーション方法、プログラム、記録媒体、及びシミュレーション装置 - Google Patents

Description

本発明は、粉粒体の時間変化をシミュレーションする方法及び装置に関する。

従来、粉粒体の時間変化のシミュレーションに、個別要素法が多く適用されていた。また、粉粒体の流動化について、粒子動力学を適用したシミュレーションが行われている。下記の非特許文献１に、薄い粉体層に振動を与えたときの粉体層の挙動を粒子動力学シミュレーションによって求める方法が示されている。

特許第４６６６３５７号公報

Physical Review Letters Vol. 77 (1996) pp.4166-4169, Keiko M. Aoki et al., "Spontaneous Wave Pattern Formation in Vibrated Granular Materials"

シミュレーション対象の粉粒体の粒子数が多くなると、巨視的な挙動をシミュレーションする際の計算時間が膨大になってしまう。上記特許文献１に、ハミルトニアンにくりこみ群変換を行うことで、ハミルトニアンの形を不変に保ち、粒子数を圧縮することができるくりこみ群分子動力学法が開示されている。くりこみ群変換に際し、ハミルトニアンの形が不変であるということは、繰り込まれた粉粒体の分解能は低下するが、繰り込まれた粉粒体の挙動は、元の粉粒体の挙動と相似になることを意味する。

上記特許文献１に示された分子動力学シミュレーションに適用される粒子間に働く力は、粒子間の相互作用ポテンシャルに依存する。これに対し、粒子動力学シミュレーションにおいては、粒子間に働く力は、粒子間相互作用ポテンシャルに依存する項と、エネルギの散逸を表す項とで表される。粒子間に働く力が、エネルギの散逸を表す項を含むため、粒子動力学シミュレーションに、上記特許文献１のくりこみ群変換を適用したのでは、粉粒体の挙動の相似性を保つことはできない。

本発明の目的は、粒子動力学シミュレーションの計算時間を短縮することが可能な粉粒体のシミュレーション方法を提供することである。本発明の他の目的は、粒子動力学シミュレーションの計算時間を短縮することが可能な粉粒体のシミュレーション装置を提供することである。

本発明の一観点によると、
各粒子に働く力が、粒子間の相互作用ポテンシャルφに依存するポテンシャル依存項と、エネルギの散逸に依存する散逸項とで表されるシミュレーション対象の粉粒体Ｓを定義し、
前記粉粒体Ｓの各粒子の運動エネルギと、前記相互作用ポテンシャルφに基づくポテンシャルエネルギとで表されるハミルトニアンの形が不変となるように、前記粉粒体Ｓの各粒子の運動方程式の前記ポテンシャル依存項に含まれる物理量をくりこみ変換し、
前記物理量のくりこみ変換による前記ポテンシャル依存項の変化率が、前記散逸項の変化率と等しくなるように、前記散逸項に含まれる物理量をくりこみ変換し、
くりこまれた粉粒体Ｓ’の各粒子の位置ベクトル及び運動量ベクトルに初期値を与えて、ある時間刻み幅で、くりこまれた前記物理量に基づいて、個々の粒子の運動方程式を数値積分することにより、くりこまれた粉粒体Ｓ’の時間発展を求める粉粒体のシミュレーション方法であって、
前記相互作用ポテンシャルφが、粒子間距離に依存する無次元関数ｆと、エネルギの次元をもつ相互作用係数εとの積εｆで表され、前記散逸項が、減衰係数ｃと、
２つの粒子の相対速度との積で表され、前記粉粒体Ｓの粒子数をＮ、各粒子の質量をｍで表し、前記粉粒体Ｓが含まれる空間の次元数をｄで表し、くりこみ回数に依存するくりこみ因子をαで表したとき、
繰り込まれた粉粒体Ｓ’の粒子数Ｎ’、各粒子の質量ｍ’、減衰係数ｃ’が、

となるように、くりこみ変換を行う粉粒体のシミュレーション方法が提供される。

本発明の他の観点によると、
各粒子に働く力が、粒子間の相互作用ポテンシャルφに依存するポテンシャル依存項と、エネルギの散逸に依存する散逸項とで表されるシミュレーション対象の粉粒体Ｓの時間発展をシミュレーションする装置であって、
入力装置と、
出力装置と、
処理ユニットと
を有し、
前記処理ユニットは、
前記入力装置を通して入力された前記粉粒体Ｓを定義する物理量を、前記粉粒体Ｓの各粒子の運動エネルギと、相互作用ポテンシャルφに基づくポテンシャルエネルギとで表されるハミルトニアンの形が不変となるように、前記粉粒体Ｓの各粒子の運動方程式のポテンシャル依存項に含まれる前記物理量をくりこみ変換し、
前記物理量のくりこみ変換による前記ポテンシャル依存項の変化率が、前記散逸項の変化率と等しくなるように、前記散逸項に含まれる物理量をくりこみ変換し、
くりこまれた前記物理量に基づいて、個々の粒子の運動方程式を、ある時間刻み幅で数値積分することにより、くりこまれた粉粒体Ｓ’の時間発展を求め、
求められた粉粒体Ｓ’の時間発展を、前記出力装置に出力するシミュレーション装置であって、
前記相互作用ポテンシャルφが、粒子間距離に依存する無次元関数ｆ、エネルギの次元をもつ相互作用係数ε、前記粉粒体Ｓを特徴づけるパラメータｒ _０ 、σ、及び粒子間距離ｒを用いて、

と表すことができ、
前記散逸項が、減衰係数ｃと、２つの粒子の相対速度との積で表され、前記粉粒体Ｓの粒子数をＮ、各粒子の質量をｍで表し、前記粉粒体Ｓが含まれる空間の次元数をｄで表し、くりこみ回数に依存するくりこみ因子をαで表したとき、
繰り込まれた粉粒体Ｓ’の粒子数Ｎ’、各粒子の質量ｍ’、相互作用係数ε’、粉粒体Ｓ’を特徴づけるパラメータσ’、及び減衰係数ｃ’が、

となるように、くりこみ変換を行い、くりこまれた粉粒体Ｓ’の粒子間の相互作用ポテンシャルφ’を、

として、くりこまれた粉粒体Ｓ’の個々の粒子の運動方程式を数値積分するシミュレーション装置が提供される。

物理量のくりこみ変換によるポテンシャル依存項の変化率が、散逸項の変化率と等しくなるように、散逸項に含まれる物理量のくりこみ変換を行うため、散逸項を含む粉粒体の時間発展をシミュレーションすることができる。このとき、粉粒体Ｓの各粒子の運動エネルギと、粉粒体Ｓの相互作用ポテンシャルφに基づくポテンシャルエネルギとで表されるハミルトニアンの形が、くりこみ変換によって不変であるため、シミュレーション対象の粉粒体Ｓの時間発展と、くりこまれた粉粒体Ｓ’の時間発展とが、相似の関係を持つ。

図１は、実施例によるシミュレーション方法のフローチャートである。 図２は、実施例によるシミュレーション方法を実行するシミュレーション装置（コンピュータ）、及びこのシミュレーション装置にシミュレーション方法を実行させるプログラムを記録した記録媒体のブロック図である。 図３Ａは、シミュレーション対象の粉粒体Ｓの初期状態を示す模式図であり、図３Ｂは、くりこみ変換を行わない粉粒体Ｓの時間展開をシミュレーションした結果を示す模式図であり、図３Ｃは、くりこまれた粉粒体Ｓ’の時間展開をシミュレーションした結果を示す模式図である。 図４は、粒子ｉ、粒子ｊ、及び粒子ｋがこの順番に一次元状に並んで配置されている状態を示す模式図である。 図５は、相互作用ポテンシャルφがレナードジョーンズ型ポテンシャルである場合の式（２１）の算出結果、及び式（２２）の数値積分結果を示すグラフである。 図６は、相互作用ポテンシャルφがモース型ポテンシャルである場合の式（２１）の算出結果、及び式（２２）の数値積分結果を示すグラフである。 図７Ａ、図７Ｂ、及び図７Ｃは、二次元格子を一次元格子に帰着させるポテンシャルムービングの方法について説明する模式図である。

実施例によるシミュレーション方法を説明する前に、実施例によるシミュレーション方法で適用される粒子動力学について簡単に説明する。粒子数Ｎのシミュレーション対象の粉粒体Ｓの粒子ｊに働く力Ｆ_ｊは、以下の式で表される。

ここで、φは、粒子間の相互作用ポテンシャルを表し、ベクトルｑ_ｊは、粒子ｊの位置ベクトル（座標）を表し、ベクトルｖ_ｊは、粒子ｊの速度ベクトルを表し、ｃは、減衰係数を表す。Σは、粒子ｊと接触している粒子について合計することを意味する。式（１）の右辺の大かっこ内の第１項は、粒子間の相互作用ポテンシャルφに依存するポテンシャル依存項であり、第２項は、エネルギの散逸に依存する散逸項である。散逸項は、減衰係数ｃと、相互作用する２つの粒子の相対速度との積で表される。

上述の式（１）から粒子ｊの運動方程式が得られる。個々の粒子に対して、運動方程式を数値積分することにより、個々の粒子の位置ベクトルｑ_ｊ、及び運動量ベクトルｐ_ｊの時間変化が求まる。

相互作用ポテンシャルφは、粒子間距離ｒに依存する無次元関数ｆと、エネルギの次元をもつ相互作用係数εとの積εｆで表される。相互作用ポテンシャルφは、一例として、以下の式で与えられる。

ここで、ｒは粒子間の距離を表し、σは粉粒体Ｓを特徴づけるパラメータである。式（２）のｒ＜ｒ_０のときの相互作用ポテンシャルφの大かっこ内の第１項及び第２項は、レナードジョーンズ型ポテンシャルに相当する。第３項の定数１／４は、ｒ＝ｒ_０において、相互作用ポテンシャルφを滑らかに接続するためのものである。

次に、図１〜図３Ｃを参照して、実施例によるシミュレーション方法について説明する。
図１に、本実施例によるシミュレーション方法のフローチャートを示す。ステップ１０において、シミュレーション対象の粉粒体Ｓを定義する。具体的には、粉粒体Ｓを特徴づける物理量である粒子数Ｎ、各粒子の質量ｍ、相互作用ポテンシャルφ、減衰係数ｃを決定する。

ステップ１１において、粉粒体Ｓの各粒子の運動エネルギと、相互作用ポテンシャルφに基づくポテンシャルエネルギとで表されるハミルトニアンの形が不変となるように、粉粒体Ｓの各粒子の運動方程式のポテンシャル依存項に含まれる物理量をくりこみ変換する。このくりこみ変換則については、特許第４６６６３５７号公報に説明されている。ハミルトニアンの形が不変であるということは、くりこみ前の粉粒体と、くりこまれた粉粒体との運動が、相互に相似であることを意味する。

一例として、くりこまれた粉粒体Ｓ’の粒子数Ｎ’、各粒子の質量ｍ’、各粒子の位置ベクトルｑ’、運動量ベクトルｐ’は、以下のくりこみ変換則によって求まる。

ここで、ｄは粉粒体Ｓが含まれる空間の次元数を表し、αはくりこみ回数に依存するくりこみ因子を表す。くりこみ回数をｎで表すと、くりこみ因子αは以下の式で表される。

式（３）のくりこみ変換則によってくりこまれた粉粒体Ｓ’に適用される相互作用ポテンシャルφは、シミュレーション対象の粉粒体Ｓの相互作用ポテンシャルφと同一である。例えば、くりこみ変換の前後において、式（２）に含まれる相互作用係数ε、パラメータσは不変である。相互作用ポテンシャルφが不変であるということは、元の粉粒体Ｓの安定状態における粒子間距離と、くりこまれた粉粒体Ｓ’の安定状態における粒子間距離とが等しいことを意味する。このため、くりこまれた粉粒体Ｓ’の巨視的な寸法が、元の粉粒体Ｓの巨視的な寸法の１／α倍に縮小される。

次に、ステップ１２において、物理量のくりこみ変換によるポテンシャル依存項の変化率と、散逸項の変化率とが等しくなるように、散逸項に含まれる物理量、具体的には減衰係数ｃ及び速度ｖのくりこみ変換を行う。以下、散逸項に含まれる物理量のくりこみ変換則について説明する。

くりこみ変換によって相互作用ポテンシャルφは不変であるため、運動方程式のポテンシャル依存項も不変である。すなわち、物理量のくりこみ変換によるポテンシャル依存項の変化率は１である。式（１）に示すように、運動方程式の散逸項は、減衰係数ｃと速度ｖとの積で表される。式（３）に示すように、くりこみ変換によって、運動量ベクトルｐの大きさがα倍になり、質量ｍがα^２倍になる。速度ｖは、ｐ／ｍに等しいため、速度ｖは、くりこみ変換によって１／α倍になる。減衰係数ｃと速度ｖとの積で表される散逸項が、くりこみ変換によって不変になるためには、粉粒体Ｓの減衰係数ｃと、くりこまれた粉粒体Ｓ’の減衰係数ｃ’との間に、以下のくりこみ変換則を適用すればよい。

次に、ステップ１３において、ステップ１１及びステップ１２でくりこまれた物理量に基づいて、くりこまれた粉粒体Ｓ’の時間発展をシミュレーションする。具体的には、粉粒体Ｓ’の各粒子の位置ベクトルｑ及び運動量ベクトルｐの初期値を設定する。くりこまれた物理量に基づいて、ある時間刻み幅で個々の粒子の運動方程式を数値積分する。この数値積分により、個々の粒子の位置ベクトルｑ及び運動量ベクトルｐの時間変化が求まる。

次に、ステップ１４において、シミュレーション結果を出力する。シミュレーション結果は、例えば各時刻における粒子の位置を、画像として表示することができる。または、粉粒体Ｓ’の巨視的な外形を画像として表示することも可能である。

図２に、実施例によるシミュレーション方法を実行するシミュレーション装置（コンピュータ）、及びこのシミュレーション装置に、図１の各ステップを実行させるプログラムを記録した記録媒体のブロック図を示す。

シミュレーション装置２０は、読取装置２１、処理ユニット２２、入力装置２３、記憶装置２４、及び出力装置２５を含む。コンピュータ読み取り可能な記録媒体３０に、図１に示したシミュレーション方法をシミュレーション装置２０に実行させるためのプログラムが記録されている。

読取装置２１が、記録媒体３０に記録されているプログラムを読み取る。読み取られたプログラムは、例えば記憶装置２４に格納され、処理ユニット２２によって実行される。以下、図１に示した各ステップ１０〜１４における具体的な手順について説明する。

ステップ１０（図１）において、シミュレーション対象の粉粒体Ｓを定義する物理量が、入力装置２３を通して取得される。これにより、シミュレーション対象の粉粒体Ｓが定義される。

ステップ１１（図１）において、くりこみ変換の回数ｎが、入力装置２３を通して取得される。処理ユニット２２は、式（４）を用いてくりこみ因子αを算出する。算出されたくりこみ因子αを用いて、運動方程式のポテンシャル依存項に含まれる物理量のくりこみ変換を行う。このくりこみ変換は、式（３）に基づいて行うことができる。ステップ１２において、運動方程式の散逸項に含まれる物理量のくりこみ変換を行う。このくりこみ変換は、式（５）を用いて行うことができる。

ステップ１３において、処理ユニット２２が、くりこまれた粉粒体Ｓ’についてシミュレーションを実行する。ステップ１４において、処理ユニット２２が、シミュレーション結果を出力装置２５に出力する。

次に、図３Ａ〜図３Ｃを参照して、上記実施例による方法でシミュレーションを行った結果と、くりこみ変換を行わないでシミュレーションを行った結果とを比較して説明する。

図３Ａに、シミュレーション対象の粉粒体Ｓの初期状態を示す。仮想的な壁４０の中に、粉粒体Ｓが収納されている。仮想的な壁４０は、鉛直方向に延びる直方体の側面を構成する。粉粒体Ｓの各粒子の質量ｍを６．３８×１０^２１／Ｎ_Ａ［ｇ］、相互作用係数εを６．８４７４×１０^２２［Ｋ］、相互作用ポテンシャルφを決定するパラメータσを３．０×１０^−３［ｍ］、減衰係数ｃを４．２４５３×１０^４［ｋｇ／ｓ］とした。ここで、Ｎ_Ａはアボガドロ定数である。空間の次元数ｄを３とした。

シミュレーション対象の粉粒体Ｓの各粒子の運動方程式に含まれる物理量に対して１回のくりこみ変換を行った。すなわち、くりこみ回数ｎが１であり、くりこみ因子αが２である。シミュレーション対象の粉粒体Ｓの粒子数が約１７万個であり、くりこまれた粉粒体Ｓ’の粒子数が約２．１万個である。シミュレーションの開始時刻において、仮想的な壁４０を取り外し、その後の粉粒体Ｓ及びＳ’の時間展開をシミュレーションにより求めた。粉粒体Ｓ及びＳ’が定常状態になった時点でシミュレーションを終了した。

図３Ｂに、くりこみ変換を行わない粉粒体Ｓの時間展開をシミュレーションした結果を示し、図３Ｃに、くりこまれた粉粒体Ｓ’の時間展開をシミュレーションした結果を示す。くりこまれた粉粒体Ｓ’の巨視的な寸法は、元の粉粒体Ｓの巨視的な寸法の１／２に縮小され、体積は１／８に縮小されている。図３Ｂと図３Ｃとを比較すると、両者の巨視的な形状はほぼ相似であることがわかる。従って、くりこまれた粉粒体Ｓ’の時間展開のシミュレーション結果から、元の粉粒体Ｓの時間展開を推測することができる。

図３Ｂに示した粉粒体Ｓのシミュレーションに要した計算時間は約１２１時間であった。これに対し、図３Ｃに示したくりこまれた粉粒体Ｓ’のシミュレーションに要した計算時間は約２６時間であった。このように、シミュレーション対象の粉粒体Ｓに対してくりこみ変換を行うことにより、シミュレーションに要する計算時間を短縮することができる。くりこみ回数を増やすと、計算時間をより短縮することが可能である。

次に、図４〜図７Ｃを参照して、他の実施例によるシミュレーション方法について説明する。以下、図１〜図３Ｃに示した実施例との相違点について説明し、共通の構成については説明を省略する。本実施例においても、粉粒体Ｓの粒子ｊに加わる力Ｆｊは、上述の式（１）で表される。

図１〜図３Ｃに示した実施例では、元の粉粒体Ｓとくりこまれた粉粒体Ｓ’とで、相互作用ポテンシャルφが同一であり、くりこまれた粉粒体Ｓ’の巨視的な寸法が元の粉粒体Ｓの巨視的な寸法の１／α倍に縮小された。以下に説明する実施例では、くりこまれた粉粒体Ｓ’の巨視的な寸法と元の粉粒体Ｓの巨視的な寸法とが同一であり、相互作用ポテンシャルφが異なる。ただし、本実施例においても、粉粒体Ｓの各粒子の運動エネルギと、相互作用ポテンシャルφに基づくポテンシャルエネルギとで表されるハミルトニアンの形は、くりこみ変換によって不変である。

以下、本実施例で適用されるくりこみ変換則について説明する。シミュレーション対象の粉粒体Ｓに適用される相互作用ポテンシャルφが、粒子間距離ｒに依存する無次元関数ｆと、相互作用係数εとを用いて以下の式で定義される。

上述の式（２）は、下記のように書き直すことができる。

従って、式（２）で表される相互作用ポテンシャルφも、式（６）で示す相互作用ポテンシャルφに含まれる。このため、相互作用ポテンシャルが式（２）で表される粉粒体Ｓも、本実施例で取り扱うことができる。

ステップ１０（図１）において定義される粉粒体Ｓを特徴づける物理量として、図１〜図３Ｃに示した実施例で粉粒体Ｓが定義された物理量に加えて、式（６）のパラメータｒ_０が定義される。

ステップ１１（図１）において、下記のくりこみ変換則により、粉粒体Ｓを定義する物理量、である粒子数Ｎ、粒子の質量ｍ、相互作用係数ε、無次元関数ｆを規定するパラメータｒ_０及びσが、それぞれ粉粒体Ｓ’を定義する物理量である粒子数Ｎ’、粒子の質量ｍ’、相互作用係数ε’、無次元関数ｆを規定するパラメータｒ_０’及びσ’に変換される。

くりこまれた粉粒体Ｓ’の相互作用ポテンシャルφ’は、以下の式で表すことができる。無次元関数ｆの形は不変である。

くりこまれた粉粒体Ｓ’のハミルトニアンＨ’は、後に詳細に説明するように、以下の式で表すことができる。ハミルトニアンの形は、くりこみ変換によって不変である。

式（８）のくりこみ変換測を適用することにより、粒子の個数が１／α^ｄ倍になり、粒子の質量がα^ｄ倍になる。このため、粉粒体Ｓの全体の質量と、粉粒体Ｓ’の全体の質量とは同一である。さらに、パラメータｒ_０、すなわち安定状態における粒子間距離がα倍になる。このため、粉粒体Ｓとくりこまれた粉粒体Ｓ’との巨視的な寸法も同一である。くりこみ変換処理の前後で、粉粒体の全質量及び寸法が不変であるため、粉粒体の密度も不変である。

次に、ステップ１２（図１）で行われる散逸項に含まれる物理量のくりこみ変換について説明する。運動方程式のポテンシャル依存項は、以下のように変形される。

ここで、Ｄｆは以下の式で定義される。

くりこまれた粉粒体Ｓ’に含まれる各粒子の運動方程式のポテンシャル依存項は、以下のように変形される。

式（１１）と式（１３）とから、以下の式が導出される。

式（１４）から、くりこみ変換によって、ポテンシャル依存項がα^{（ｄ−１）}倍になることがわかる。

式（８）の質量ｍのくりこみ変換則、及び式（１０）の運動量のくりこみ変換則から、くりこみ変換によって速度ｖは不変であることがわかる。式（１）の散逸項に含まれる速度ｖが不変であるから、散逸項をα^{（ｄ−１）}倍にするためには、粉粒体Ｓの減衰係数ｃから、粉粒体Ｓ’の減衰係数ｃ’へのくりこみ変換則を以下の式で定義すればよい。

本実施例によるステップ１３及びステップ１４（図１）の処理は、図１〜図３Ｃに示した処理と共通である。

次に、図４〜図７を参照して、くりこまれた粉粒体Ｓ’のハミルトニアンＨ’（式（１０））の導出について説明する。粉粒体Ｓに対する分配関数Ｚ（β）の積分の一部を実行し、ハミルトニアンを粗視化することにより、くりこまれた粉粒体Ｓ’のハミルトニアンＨ’が得られる。

粒子数が一定の正準集団（Canonical ensemble）に対する分配関数Ｚ（β）は、以下の式で表される。

ここで、ｄΓ_Ｎは位相空間内の体積要素であり、下記の式で表される。

ここで、ｈはプランク定数である。Ｗ_Ｎは、すべての状態の量子論的な和と、位相空間に亘る積分とが一致するように決められる。

まず、粒子間の相互作用ポテンシャルの粗視化について説明し、次に運動エネルギの粗視化について説明する。続いて、相互作用ポテンシャルの粗視化、及び運動エネルギの粗視化を踏まえて、くりこみ変換測を定義する。

［粒子間の相互作用ポテンシャルの粗視化］
まず、一次元鎖状に配置された粉粒体における相互作用ポテンシャルの粗視化について説明する。その後、単純立方格子状に配置された粉粒体における相互作用ポテンシャルについて説明する。

図４に示すように、粒子ｉ、粒子ｊ、及び粒子ｋがこの順番に一次元状に並んで配置されている。粒子ｊが関与する相互作用を書き出し、粒子ｉと粒子ｋとの中間にある粒子ｊの位置座標について積分を実行することにより、相互作用ポテンシャルの粗視化を行うことができる。まず、第二近接以上の粒子からの寄与を取り込むために、ポテンシャルムービング（Potential Moving）の方法を用いることができる。ポテンシャルムービングの方法については、Leo P. Kadanoff, STATISTICAL PHYSICS Static, Dynamics and Renormalization, Chap. 14, World Scientific (1999)に解説されている。

第二近接以上の粒子からの寄与を取り込んだ相互作用ポテンシャルφティルダは、下記の式で表すことができる。

ここで、ａは平衡状態における粒子間距離を表す。平衡状態における粒子間距離ａは、相互作用ポテンシャルφが最小となる距離であるｒ_０に等しいと近似することができる。

複数の粒子が一次元状に配置されているため、粒子ｊの位置ベクトルｑ_ｊは、一次元座標で表すことができる。粒子ｊの位置をｑ_ｊで表すと、粒子ｊに対して、最近接の粒子によって作られるケージポテンシャルは、下記の式で表される。

位置ベクトルｑ_ｊについて積分を実行すれば、式（１９）を用いて下記の式が得られる。

ここで、ｒ_ａは粒子の直径を表し、ｚ（ｑ_ｉ−ｑ_ｋ）及びＰ（ｑ_ｉ−ｑ_ｋ）は、下記の式で表される。

積分領域はケージポテンシャルの内部領域に制限される。

次に、ｚ（ｑ_ｉ−ｑ_ｋ）を具体的に計算する。相互作用ポテンシャルφがレナードジョーンズ型ポテンシャルの場合、φ^（２ｎ）は、下記の式で表される。

相互作用ポテンシャルφがモース型ポテンシャルの場合、φ^（２ｎ）は、下記の式で表される。

式（２３）または式（２４）を式（２２）に代入して数値積分を行う。式（２２）に式（２３）または式（２４）を代入する際には、式（１８）を用いる。数値積分において、式（２２）の積分範囲に現れているａは、近似的にｒ_０に等しいと仮定する。

図５に、相互作用ポテンシャルφがレナードジョーンズ型ポテンシャルである場合の式（２１）の算出結果、及び式（２２）の数値積分結果を示す。図６に、相互作用ポテンシャルφがモース型ポテンシャルである場合の式（２１）の算出結果、及び式（２２）の数値積分結果を示す。粒子ｉ、粒子ｊ、粒子ｋ、粒子ｌ、及び粒子ｍがこの順番に一次元状に配置されている場合の粒子ｋの位置座標がｑ_ｋで表される。図５及び図６の横軸は、ｑ_ｋ／２を表し、縦軸は、Ｐ（ｑ_ｉ−ｑ_ｋ）Ｐ（ｑ_ｋ−ｑ_ｍ）、及びｚ（ｑ_ｉ−ｑ_ｋ）ｚ（ｑ_ｋ−ｑ_ｍ）を、対数目盛で表す。図５に至る数値計算において、ε／ｋ_ＢＴ＝２．０、ｒ_０／σ＝１．１２とした。図６に至る数値計算において、ε／ｋ_ＢＴ＝２．０、ｒ_０／σ＝２．２４とした。

相互作用ポテンシャルφがレナードジョーンズ型ポテンシャルである場合、及びモース型ポテンシャルである場合のいずれにおいても、Ｐ（ｑ_ｉ−ｑ_ｋ）Ｐ（ｑ_ｋ−ｑ_ｍ）の変化に比べて、ｚ（ｑ_ｉ−ｑ_ｋ）ｚ（ｑ_ｋ−ｑ_ｍ）の変化が緩やかであることがわかる。このため、Ｐ（ｑ_ｉ−ｑ_ｋ）Ｐ（ｑ_ｋ−ｑ_ｍ）に対してｚ（ｑ_ｉ−ｑ_ｋ）ｚ（ｑ_ｋ−ｑ_ｍ）は、ほぼ定数であると近似できる。

粒子ｋが位置座標ｑ_ｋに存在する確率ｐ（ｑ_ｋ）は、以下のように近似することができる。

従って、下記の式が導出される。

ここまでは、一次元状に複数の粒子が配置された粉粒体の相互作用ポテンシャルの粗視化について説明した。多次元の粉粒体の相互作用ポテンシャルは、ポテンシャルムービングの方法により実現できる。

図７Ａ〜図７Ｃを参照して、二次元格子を一次元格子に帰着させるポテンシャルムービングの方法について説明する。

図７Ａに示すように、二次元正方格子の格子点の位置に粒子が配置されている。最近接する粒子間の相互作用（ニアレストネイバーカップリング）を実線で示す。粒子が配列する１つの方向をｘ方向とし、それに直交する方向をｙ方向とする。

図７Ｂに示すように、ｘ方向に並んだ粒子のうち、１つ置きに配置された粒子（中空の円で表された粒子）の変位について積分を実行することを考える。積分の対象となる粒子（消去したい粒子）を被積分粒子と呼ぶこととする。

図７Ｃに示すように、ポテンシャルムービングの方法では、被積分粒子同士のニアレストネイバーカップリングを、ｘ方向に関して両隣に配置された粒子に振り分ける。振り分けられた相互作用が加算された粒子間相互作用（ダブルカップリング）を二重線で表す。二重線で表されたダブルカップリングは、元のニアレストネイバーカップリングの２倍の強さになる。この手法を用いて、二次元格子を一次元鎖に変換することができる。一次元鎖の粉粒体においては、図４を参照して説明した手法により相互作用ポテンシャルの粗視化を行うことができる。シミュレーション対象の粉粒体が三次元格子を構成する場合には、二次元格子を一次元鎖に変換する手順をｘ方向、ｙ方向、及びｚ方向の３方向に関して繰り返せばよい。このようにして、多次元格子を構成する粉粒体の粗視化を行うことができる。

次元数ｄの多次元格子を構成する粉粒体の相互作用ポテンシャルの粗視化は、下記の式で表される。

ここで、＜ｉ，ｊ＞は、最近接格子間で和を取ることを意味する。

すべての相互作用の和に変更すると、下記の式が得られる。

［運動エネルギの粗視化］
次に、運動エネルギの粗視化について説明する。運動エネルギについては、容易に積分を実行することができ、下記の式が導出される。

式（２９）の導出に際し、下記の式が利用される。ここで、運動量ベクトルｐ_ｊ ^２は、ベクトルの内積を意味する。

［くりこみ変換測の導出］
次に、上述の相互作用ポテンシャルの粗視化、及び運動エネルギの粗視化から導出されるくりこみ変換測について説明する。

式（１６）に、式（２８）及び式（２９）を代入することにより、結果に影響しない係数を除いて下記の式が得られる。

式（３１）から、粗視化されたハミルトニアン（くりこまれた粉粒体Ｓ’のハミルトニアン）Ｈ’は、以下の式で表される。

ハミルトニアンを粗視化する際の結合定数のリストをＫで表す。結合定数のリストＫは下記のように表される。

くりこみ変換Ｒを以下のように定義する。

ｎ回のくりこみ変換実行後の結合係数のリストＫ_ｎは、下記の式で表される。

従って、ｎ回のくりこみ変換が行われたハミルトニアンＨｎは、下記の式で表される。

ここで、くりこまれた粉粒体Ｓ’における運動量ベクトルｐ_ｊ’は、下記の式で表される。

式（３６）から、式（８）に示したくりこみ変換測、及び式（１０）に示したくりこまれた粉粒体Ｓ’のハミルトニアンＨ’が導き出される。

本実施例においても、くりこまれた粉粒体Ｓ’の各粒子の運動エネルギと、相互作用ポテンシャルφに基づくポテンシャルエネルギとで表されるハミルトニアンの形が、くりこみ前の粉粒体Ｓのハミルトニアンの形と同一である。このため、本実施例においても、図１〜図３Ｃに示した実施例と同様に、くりこみ前の粉粒体Ｓと、くりこまれた粉粒体Ｓ’との運動が、相互に相似である。

以上実施例に沿って本発明を説明したが、本発明はこれらに制限されるものではない。例えば、種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。

１０、１１、１２、１３、１４ ステップ
２０ シミュレーション装置（コンピュータ）
２１ 読取装置
２２ 処理ユニット
２３ 入力装置
２４ 記憶装置
２５ 出力装置
３０ 記録媒体
４０ 仮想的な壁
Ｈ’、Ｈｎ ハミルトニアン
Ｎ、Ｎ’ 粒子数
Ｓ シミュレーション対象の粉粒体
Ｓ’ 繰り込まれた粉粒体
ｃ、ｃ’ 減衰係数
ｄ 次元数
ｆ 無次元関数
ｍ、ｍ’ 質量
ｐ、ｐ’ 運動量ベクトル
ｑ、ｑ’ 位置ベクトル
ｒ 粒子間距離
ｖ 速度
α くりこみ因子
ε、ε’ 相互作用係数
φ、φ’ 相互作用ポテンシャル

Claims (8)

1. 各粒子に働く力が、粒子間の相互作用ポテンシャルφに依存するポテンシャル依存項と、エネルギの散逸に依存する散逸項とで表されるシミュレーション対象の粉粒体Ｓを定義し、
   前記粉粒体Ｓの各粒子の運動エネルギと、前記相互作用ポテンシャルφに基づくポテンシャルエネルギとで表されるハミルトニアンの形が不変となるように、前記粉粒体Ｓの各粒子の運動方程式の前記ポテンシャル依存項に含まれる物理量をくりこみ変換し、
   前記物理量のくりこみ変換による前記ポテンシャル依存項の変化率が、前記散逸項の変化率と等しくなるように、前記散逸項に含まれる物理量をくりこみ変換し、
   くりこまれた粉粒体Ｓ’の各粒子の位置ベクトル及び運動量ベクトルに初期値を与えて、ある時間刻み幅で、くりこまれた前記物理量に基づいて、個々の粒子の運動方程式を数値積分することにより、くりこまれた粉粒体Ｓ’の時間発展を求める粉粒体のシミュレーション方法であって、
   前記相互作用ポテンシャルφが、粒子間距離に依存する無次元関数ｆと、エネルギの次元をもつ相互作用係数εとの積εｆで表され、前記散逸項が、減衰係数ｃと、
   ２つの粒子の相対速度との積で表され、前記粉粒体Ｓの粒子数をＮ、各粒子の質量をｍで表し、前記粉粒体Ｓが含まれる空間の次元数をｄで表し、くりこみ回数に依存するくりこみ因子をαで表したとき、
   繰り込まれた粉粒体Ｓ’の粒子数Ｎ’、各粒子の質量ｍ’、減衰係数ｃ’が、
   

   

   となるように、くりこみ変換を行う粉粒体のシミュレーション方法。
2. 各粒子に働く力が、粒子間の相互作用ポテンシャルφに依存するポテンシャル依存項と、エネルギの散逸に依存する散逸項とで表されるシミュレーション対象の粉粒体Ｓを定義し、
   前記粉粒体Ｓの各粒子の運動エネルギと、前記相互作用ポテンシャルφに基づくポテンシャルエネルギとで表されるハミルトニアンの形が不変となるように、前記粉粒体Ｓの各粒子の運動方程式の前記ポテンシャル依存項に含まれる物理量をくりこみ変換し、
   前記物理量のくりこみ変換による前記ポテンシャル依存項の変化率が、前記散逸項の変化率と等しくなるように、前記散逸項に含まれる物理量をくりこみ変換し、
   くりこまれた粉粒体Ｓ’の各粒子の位置ベクトル及び運動量ベクトルに初期値を与えて、ある時間刻み幅で、くりこまれた前記物理量に基づいて、個々の粒子の運動方程式を数値積分することにより、くりこまれた粉粒体Ｓ’の時間発展を求める粉粒体のシミュレーション方法であって、
   前記相互作用ポテンシャルφが、粒子間距離に依存する無次元関数ｆ、エネルギの次元をもつ相互作用係数ε、前記粉粒体Ｓを特徴づけるパラメータｒ_０、σ、及び粒子間距離ｒを用いて、
   

   

   と表すことができ、
   前記散逸項が、減衰係数ｃと、２つの粒子の相対速度との積で表され、前記粉粒体Ｓの粒子数をＮ、各粒子の質量をｍで表し、前記粉粒体Ｓが含まれる空間の次元数をｄで表し、くりこみ回数に依存するくりこみ因子をαで表したとき、
   繰り込まれた粉粒体Ｓ’の粒子数Ｎ’、各粒子の質量ｍ’、相互作用係数ε’、粉粒体Ｓ’を特徴づけるパラメータσ’、及び減衰係数ｃ’が、
   

   

   となるように、くりこみ変換を行い、くりこまれた粉粒体Ｓ’の粒子間の相互作用ポテンシャルφ’を、
   

   

   として、くりこまれた粉粒体Ｓ’の個々の粒子の運動方程式を数値積分する請求項１に記載の粉粒体のシミュレーション方法。
3. 各粒子に働く力が、粒子間の相互作用ポテンシャルφに依存するポテンシャル依存項と、エネルギの散逸に依存する散逸項とで表されるシミュレーション対象の粉粒体Ｓを定義する物理量を取得する手順と、
   前記粉粒体Ｓの各粒子の運動エネルギと、前記相互作用ポテンシャルφに基づくポテンシャルエネルギとで表されるハミルトニアンの形が不変となるように、前記粉粒体Ｓの各粒子の運動方程式の前記ポテンシャル依存項に含まれる前記物理量をくりこみ変換する手順と、
   前記物理量のくりこみ変換による前記ポテンシャル依存項の変化率が、前記散逸項の変化率と等しくなるように、前記散逸項に含まれる物理量をくりこみ変換する手順と、
   くりこまれた粉粒体Ｓ’の各粒子の位置ベクトル及び運動量ベクトルに初期値を与えて、ある時間刻み幅で、くりこまれた前記物理量に基づいて、個々の粒子の運動方程式を数値積分することにより、くりこまれた粉粒体Ｓ’の時間発展を求める手順と
   をコンピュータに実行させるプログラムであって、
   前記相互作用ポテンシャルφが、粒子間距離に依存する無次元関数ｆと、エネルギの次元をもつ相互作用係数εとの積εｆで表され、前記散逸項が、減衰係数ｃと、
   ２つの粒子の相対速度との積で表され、前記粉粒体Ｓの粒子数をＮ、各粒子の質量をｍで表し、前記粉粒体Ｓが含まれる空間の次元数をｄで表し、くりこみ回数に依存するくりこみ因子をαで表したとき、
   繰り込まれた粉粒体Ｓ’の粒子数Ｎ’、各粒子の質量ｍ’、減衰係数ｃ’が、
   

   

   となるように、くりこみ変換を行うプログラム。
4. 各粒子に働く力が、粒子間の相互作用ポテンシャルφに依存するポテンシャル依存項と、エネルギの散逸に依存する散逸項とで表されるシミュレーション対象の粉粒体Ｓを定義する物理量を取得する手順と、
   前記粉粒体Ｓの各粒子の運動エネルギと、前記相互作用ポテンシャルφに基づくポテンシャルエネルギとで表されるハミルトニアンの形が不変となるように、前記粉粒体Ｓの各粒子の運動方程式の前記ポテンシャル依存項に含まれる前記物理量をくりこみ変換する手順と、
   前記物理量のくりこみ変換による前記ポテンシャル依存項の変化率が、前記散逸項の変化率と等しくなるように、前記散逸項に含まれる物理量をくりこみ変換する手順と、
   くりこまれた粉粒体Ｓ’の各粒子の位置ベクトル及び運動量ベクトルに初期値を与えて、ある時間刻み幅で、くりこまれた前記物理量に基づいて、個々の粒子の運動方程式を数値積分することにより、くりこまれた粉粒体Ｓ’の時間発展を求める手順と
   をコンピュータに実行させるプログラムであって、
   前記相互作用ポテンシャルφが、粒子間距離に依存する無次元関数ｆ、エネルギの次元をもつ相互作用係数ε、前記粉粒体Ｓを特徴づけるパラメータｒ _０ 、σ、及び粒子間距離ｒを用いて、
   

   

   と表すことができ、
   前記散逸項が、減衰係数ｃと、２つの粒子の相対速度との積で表され、前記粉粒体Ｓの粒子数をＮ、各粒子の質量をｍで表し、前記粉粒体Ｓが含まれる空間の次元数をｄで表し、くりこみ回数に依存するくりこみ因子をαで表したとき、
   繰り込まれた粉粒体Ｓ’の粒子数Ｎ’、各粒子の質量ｍ’、相互作用係数ε’、粉粒体Ｓ’を特徴づけるパラメータσ’、及び減衰係数ｃ’が、
   

   

   となるように、くりこみ変換を行い、くりこまれた粉粒体Ｓ’の粒子間の相互作用ポテンシャルφ’を、
   

   

   として、くりこまれた粉粒体Ｓ’の個々の粒子の運動方程式を数値積分するプログラム。
5. 各粒子に働く力が、粒子間の相互作用ポテンシャルφに依存するポテンシャル依存項と、エネルギの散逸に依存する散逸項とで表されるシミュレーション対象の粉粒体Ｓを定義する物理量を取得する手順と、
   前記粉粒体Ｓの各粒子の運動エネルギと、前記相互作用ポテンシャルφに基づくポテンシャルエネルギとで表されるハミルトニアンの形が不変となるように、前記粉粒体Ｓの各粒子の運動方程式の前記ポテンシャル依存項に含まれる前記物理量をくりこみ変換する手順と、
   前記物理量のくりこみ変換による前記ポテンシャル依存項の変化率が、前記散逸項の変化率と等しくなるように、前記散逸項に含まれる物理量をくりこみ変換する手順と、
   くりこまれた粉粒体Ｓ’の各粒子の位置ベクトル及び運動量ベクトルに初期値を与えて、ある時間刻み幅で、くりこまれた前記物理量に基づいて、個々の粒子の運動方程式を数値積分することにより、くりこまれた粉粒体Ｓ’の時間発展を求める手順と
   をコンピュータに実行させるプログラムを記録したコンピュータ読み取り可能な記録媒体であって、
   前記相互作用ポテンシャルφが、粒子間距離に依存する無次元関数ｆと、エネルギの次元をもつ相互作用係数εとの積εｆで表され、前記散逸項が、減衰係数ｃと、
   ２つの粒子の相対速度との積で表され、前記粉粒体Ｓの粒子数をＮ、各粒子の質量をｍで表し、前記粉粒体Ｓが含まれる空間の次元数をｄで表し、くりこみ回数に依存するくりこみ因子をαで表したとき、
   繰り込まれた粉粒体Ｓ’の粒子数Ｎ’、各粒子の質量ｍ’、減衰係数ｃ’が、
   

   

   となるように、くりこみ変換を行うプログラムを記録したコンピュータ読み取り可能な記録媒体。
6. 各粒子に働く力が、粒子間の相互作用ポテンシャルφに依存するポテンシャル依存項と、エネルギの散逸に依存する散逸項とで表されるシミュレーション対象の粉粒体Ｓを定義する物理量を取得する手順と、
   前記粉粒体Ｓの各粒子の運動エネルギと、前記相互作用ポテンシャルφに基づくポテンシャルエネルギとで表されるハミルトニアンの形が不変となるように、前記粉粒体Ｓの各粒子の運動方程式の前記ポテンシャル依存項に含まれる前記物理量をくりこみ変換する手順と、
   前記物理量のくりこみ変換による前記ポテンシャル依存項の変化率が、前記散逸項の変化率と等しくなるように、前記散逸項に含まれる物理量をくりこみ変換する手順と、
   くりこまれた粉粒体Ｓ’の各粒子の位置ベクトル及び運動量ベクトルに初期値を与えて、ある時間刻み幅で、くりこまれた前記物理量に基づいて、個々の粒子の運動方程式を数値積分することにより、くりこまれた粉粒体Ｓ’の時間発展を求める手順と
   をコンピュータに実行させるプログラムを記録したコンピュータ読み取り可能な記録媒体であって、
   前記相互作用ポテンシャルφが、粒子間距離に依存する無次元関数ｆ、エネルギの次元をもつ相互作用係数ε、前記粉粒体Ｓを特徴づけるパラメータｒ _０ 、σ、及び粒子間距離ｒを用いて、
   

   

   と表すことができ、
   前記散逸項が、減衰係数ｃと、２つの粒子の相対速度との積で表され、前記粉粒体Ｓの粒子数をＮ、各粒子の質量をｍで表し、前記粉粒体Ｓが含まれる空間の次元数をｄで表し、くりこみ回数に依存するくりこみ因子をαで表したとき、
   繰り込まれた粉粒体Ｓ’の粒子数Ｎ’、各粒子の質量ｍ’、相互作用係数ε’、粉粒体Ｓ’を特徴づけるパラメータσ’、及び減衰係数ｃ’が、
   

   

   となるように、くりこみ変換を行い、くりこまれた粉粒体Ｓ’の粒子間の相互作用ポテンシャルφ’を、
   

   

   として、くりこまれた粉粒体Ｓ’の個々の粒子の運動方程式を数値積分するプログラムを記録したコンピュータ読み取り可能な記録媒体。
7. 各粒子に働く力が、粒子間の相互作用ポテンシャルφに依存するポテンシャル依存項と、エネルギの散逸に依存する散逸項とで表されるシミュレーション対象の粉粒体Ｓの時間発展をシミュレーションする装置であって、
   入力装置と、
   出力装置と、
   処理ユニットと
   を有し、
   前記処理ユニットは、
   前記入力装置を通して入力された前記粉粒体Ｓを定義する物理量を、前記粉粒体Ｓの各粒子の運動エネルギと、相互作用ポテンシャルφに基づくポテンシャルエネルギとで表されるハミルトニアンの形が不変となるように、前記粉粒体Ｓの各粒子の運動方程式のポテンシャル依存項に含まれる前記物理量をくりこみ変換し、
   前記物理量のくりこみ変換による前記ポテンシャル依存項の変化率が、前記散逸項の変化率と等しくなるように、前記散逸項に含まれる物理量をくりこみ変換し、
   くりこまれた前記物理量に基づいて、個々の粒子の運動方程式を、ある時間刻み幅で数値積分することにより、くりこまれた粉粒体Ｓ’の時間発展を求め、
   求められた粉粒体Ｓ’の時間発展を、前記出力装置に出力するシミュレーション装置であって、
   前記相互作用ポテンシャルφが、粒子間距離に依存する無次元関数ｆと、エネルギの次元をもつ相互作用係数εとの積εｆで表され、前記散逸項が、減衰係数ｃと、
   ２つの粒子の相対速度との積で表され、前記粉粒体Ｓの粒子数をＮ、各粒子の質量をｍで表し、前記粉粒体Ｓが含まれる空間の次元数をｄで表し、くりこみ回数に依存するくりこみ因子をαで表したとき、
   繰り込まれた粉粒体Ｓ’の粒子数Ｎ’、各粒子の質量ｍ’、減衰係数ｃ’が、
   

   

   となるように、くりこみ変換を行うシミュレーション装置。
8. 各粒子に働く力が、粒子間の相互作用ポテンシャルφに依存するポテンシャル依存項と、エネルギの散逸に依存する散逸項とで表されるシミュレーション対象の粉粒体Ｓの時間発展をシミュレーションする装置であって、
   入力装置と、
   出力装置と、
   処理ユニットと
   を有し、
   前記処理ユニットは、
   前記入力装置を通して入力された前記粉粒体Ｓを定義する物理量を、前記粉粒体Ｓの各粒子の運動エネルギと、相互作用ポテンシャルφに基づくポテンシャルエネルギとで表されるハミルトニアンの形が不変となるように、前記粉粒体Ｓの各粒子の運動方程式のポテンシャル依存項に含まれる前記物理量をくりこみ変換し、
   前記物理量のくりこみ変換による前記ポテンシャル依存項の変化率が、前記散逸項の変化率と等しくなるように、前記散逸項に含まれる物理量をくりこみ変換し、
   くりこまれた前記物理量に基づいて、個々の粒子の運動方程式を、ある時間刻み幅で数値積分することにより、くりこまれた粉粒体Ｓ’の時間発展を求め、
   求められた粉粒体Ｓ’の時間発展を、前記出力装置に出力するシミュレーション装置であって、
   前記相互作用ポテンシャルφが、粒子間距離に依存する無次元関数ｆ、エネルギの次元をもつ相互作用係数ε、前記粉粒体Ｓを特徴づけるパラメータｒ _０ 、σ、及び粒子間距離ｒを用いて、
   

   

   と表すことができ、
   前記散逸項が、減衰係数ｃと、２つの粒子の相対速度との積で表され、前記粉粒体Ｓの粒子数をＮ、各粒子の質量をｍで表し、前記粉粒体Ｓが含まれる空間の次元数をｄで表し、くりこみ回数に依存するくりこみ因子をαで表したとき、
   繰り込まれた粉粒体Ｓ’の粒子数Ｎ’、各粒子の質量ｍ’、相互作用係数ε’、粉粒体Ｓ’を特徴づけるパラメータσ’、及び減衰係数ｃ’が、
   

   

   となるように、くりこみ変換を行い、くりこまれた粉粒体Ｓ’の粒子間の相互作用ポテンシャルφ’を、
   

   

   として、くりこまれた粉粒体Ｓ’の個々の粒子の運動方程式を数値積分するシミュレーション装置。
   

Images