(57)【特許請求の範囲】
【請求項１】
各粒子に働く力が、粒子間の相互作用ポテンシャルφに依存するポテンシャル依存項と、エネルギの散逸に依存する散逸項とで表されるシミュレーション対象の粉粒体Ｓを定義し、
前記粉粒体Ｓの各粒子の運動エネルギと、前記相互作用ポテンシャルφに基づくポテンシャルエネルギとで表されるハミルトニアンの形が不変となるように、前記粉粒体Ｓの各粒子の運動方程式の前記ポテンシャル依存項に含まれる物理量をくりこみ変換し、
前記物理量のくりこみ変換による前記ポテンシャル依存項の変化率が、前記散逸項の変化率と等しくなるように、前記散逸項に含まれる物理量をくりこみ変換し、
くりこまれた粉粒体Ｓ’の各粒子の位置ベクトル及び運動量ベクトルに初期値を与えて、ある時間刻み幅で、くりこまれた前記物理量に基づいて、個々の粒子の運動方程式を数値積分することにより、くりこまれた粉粒体Ｓ’の時間発展を求める粉粒体のシミュレーション方法であって、
前記相互作用ポテンシャルφが、粒子間距離に依存する無次元関数ｆと、エネルギの次元をもつ相互作用係数εとの積εｆで表され、前記散逸項が、減衰係数ｃと、
２つの粒子の相対速度との積で表され、前記粉粒体Ｓの粒子数をＮ、各粒子の質量をｍで表し、前記粉粒体Ｓが含まれる空間の次元数をｄで表し、くりこみ回数に依存するくりこみ因子をαで表したとき、
繰り込まれた粉粒体Ｓ’の粒子数Ｎ’、各粒子の質量ｍ’、減衰係数ｃ’が、
【数１】
0006472346.tif 000043
となるように、くりこみ変換を行う粉粒体のシミュレーション方法。
【請求項２】
各粒子に働く力が、粒子間の相互作用ポテンシャルφに依存するポテンシャル依存項と、エネルギの散逸に依存する散逸項とで表されるシミュレーション対象の粉粒体Ｓを定義し、
前記粉粒体Ｓの各粒子の運動エネルギと、前記相互作用ポテンシャルφに基づくポテンシャルエネルギとで表されるハミルトニアンの形が不変となるように、前記粉粒体Ｓの各粒子の運動方程式の前記ポテンシャル依存項に含まれる物理量をくりこみ変換し、
前記物理量のくりこみ変換による前記ポテンシャル依存項の変化率が、前記散逸項の変化率と等しくなるように、前記散逸項に含まれる物理量をくりこみ変換し、
くりこまれた粉粒体Ｓ’の各粒子の位置ベクトル及び運動量ベクトルに初期値を与えて、ある時間刻み幅で、くりこまれた前記物理量に基づいて、個々の粒子の運動方程式を数値積分することにより、くりこまれた粉粒体Ｓ’の時間発展を求める粉粒体のシミュレーション方法であって、
前記相互作用ポテンシャルφが、粒子間距離に依存する無次元関数ｆ、エネルギの次元をもつ相互作用係数ε、前記粉粒体Ｓを特徴づけるパラメータｒ０、σ、及び粒子間距離ｒを用いて、
【数２】
0006472346.tif 000044
と表すことができ、
前記散逸項が、減衰係数ｃと、２つの粒子の相対速度との積で表され、前記粉粒体Ｓの粒子数をＮ、各粒子の質量をｍで表し、前記粉粒体Ｓが含まれる空間の次元数をｄで表し、くりこみ回数に依存するくりこみ因子をαで表したとき、
繰り込まれた粉粒体Ｓ’の粒子数Ｎ’、各粒子の質量ｍ’、相互作用係数ε’、粉粒体Ｓ’を特徴づけるパラメータσ’、及び減衰係数ｃ’が、
【数３】
0006472346.tif 000045
となるように、くりこみ変換を行い、くりこまれた粉粒体Ｓ’の粒子間の相互作用ポテンシャルφ’を、
【数４】
0006472346.tif 000046
として、くりこまれた粉粒体Ｓ’の個々の粒子の運動方程式を数値積分する請求項１に記載の粉粒体のシミュレーション方法。
【請求項３】
各粒子に働く力が、粒子間の相互作用ポテンシャルφに依存するポテンシャル依存項と、エネルギの散逸に依存する散逸項とで表されるシミュレーション対象の粉粒体Ｓを定義する物理量を取得する手順と、
前記粉粒体Ｓの各粒子の運動エネルギと、前記相互作用ポテンシャルφに基づくポテンシャルエネルギとで表されるハミルトニアンの形が不変となるように、前記粉粒体Ｓの各粒子の運動方程式の前記ポテンシャル依存項に含まれる前記物理量をくりこみ変換する手順と、
前記物理量のくりこみ変換による前記ポテンシャル依存項の変化率が、前記散逸項の変化率と等しくなるように、前記散逸項に含まれる物理量をくりこみ変換する手順と、
くりこまれた粉粒体Ｓ’の各粒子の位置ベクトル及び運動量ベクトルに初期値を与えて、ある時間刻み幅で、くりこまれた前記物理量に基づいて、個々の粒子の運動方程式を数値積分することにより、くりこまれた粉粒体Ｓ’の時間発展を求める手順と
をコンピュータに実行させるプログラムであって、
前記相互作用ポテンシャルφが、粒子間距離に依存する無次元関数ｆと、エネルギの次元をもつ相互作用係数εとの積εｆで表され、前記散逸項が、減衰係数ｃと、
２つの粒子の相対速度との積で表され、前記粉粒体Ｓの粒子数をＮ、各粒子の質量をｍで表し、前記粉粒体Ｓが含まれる空間の次元数をｄで表し、くりこみ回数に依存するくりこみ因子をαで表したとき、
繰り込まれた粉粒体Ｓ’の粒子数Ｎ’、各粒子の質量ｍ’、減衰係数ｃ’が、
【数５】
0006472346.tif 000047
となるように、くりこみ変換を行うプログラム。
【請求項４】
各粒子に働く力が、粒子間の相互作用ポテンシャルφに依存するポテンシャル依存項と、エネルギの散逸に依存する散逸項とで表されるシミュレーション対象の粉粒体Ｓを定義する物理量を取得する手順と、
前記粉粒体Ｓの各粒子の運動エネルギと、前記相互作用ポテンシャルφに基づくポテンシャルエネルギとで表されるハミルトニアンの形が不変となるように、前記粉粒体Ｓの各粒子の運動方程式の前記ポテンシャル依存項に含まれる前記物理量をくりこみ変換する手順と、
前記物理量のくりこみ変換による前記ポテンシャル依存項の変化率が、前記散逸項の変化率と等しくなるように、前記散逸項に含まれる物理量をくりこみ変換する手順と、
くりこまれた粉粒体Ｓ’の各粒子の位置ベクトル及び運動量ベクトルに初期値を与えて、ある時間刻み幅で、くりこまれた前記物理量に基づいて、個々の粒子の運動方程式を数値積分することにより、くりこまれた粉粒体Ｓ’の時間発展を求める手順と
をコンピュータに実行させるプログラムであって、
前記相互作用ポテンシャルφが、粒子間距離に依存する無次元関数ｆ、エネルギの次元をもつ相互作用係数ε、前記粉粒体Ｓを特徴づけるパラメータｒ０、σ、及び粒子間距離ｒを用いて、
【数６】
0006472346.tif 000048
と表すことができ、
前記散逸項が、減衰係数ｃと、２つの粒子の相対速度との積で表され、前記粉粒体Ｓの粒子数をＮ、各粒子の質量をｍで表し、前記粉粒体Ｓが含まれる空間の次元数をｄで表し、くりこみ回数に依存するくりこみ因子をαで表したとき、
繰り込まれた粉粒体Ｓ’の粒子数Ｎ’、各粒子の質量ｍ’、相互作用係数ε’、粉粒体Ｓ’を特徴づけるパラメータσ’、及び減衰係数ｃ’が、
【数７】
0006472346.tif 000049
となるように、くりこみ変換を行い、くりこまれた粉粒体Ｓ’の粒子間の相互作用ポテンシャルφ’を、
【数８】
0006472346.tif 000050
として、くりこまれた粉粒体Ｓ’の個々の粒子の運動方程式を数値積分するプログラム。
【請求項５】
各粒子に働く力が、粒子間の相互作用ポテンシャルφに依存するポテンシャル依存項と、エネルギの散逸に依存する散逸項とで表されるシミュレーション対象の粉粒体Ｓを定義する物理量を取得する手順と、
前記粉粒体Ｓの各粒子の運動エネルギと、前記相互作用ポテンシャルφに基づくポテンシャルエネルギとで表されるハミルトニアンの形が不変となるように、前記粉粒体Ｓの各粒子の運動方程式の前記ポテンシャル依存項に含まれる前記物理量をくりこみ変換する手順と、
前記物理量のくりこみ変換による前記ポテンシャル依存項の変化率が、前記散逸項の変化率と等しくなるように、前記散逸項に含まれる物理量をくりこみ変換する手順と、
くりこまれた粉粒体Ｓ’の各粒子の位置ベクトル及び運動量ベクトルに初期値を与えて、ある時間刻み幅で、くりこまれた前記物理量に基づいて、個々の粒子の運動方程式を数値積分することにより、くりこまれた粉粒体Ｓ’の時間発展を求める手順と
をコンピュータに実行させるプログラムを記録したコンピュータ読み取り可能な記録媒体であって、
前記相互作用ポテンシャルφが、粒子間距離に依存する無次元関数ｆと、エネルギの次元をもつ相互作用係数εとの積εｆで表され、前記散逸項が、減衰係数ｃと、
２つの粒子の相対速度との積で表され、前記粉粒体Ｓの粒子数をＮ、各粒子の質量をｍで表し、前記粉粒体Ｓが含まれる空間の次元数をｄで表し、くりこみ回数に依存するくりこみ因子をαで表したとき、
繰り込まれた粉粒体Ｓ’の粒子数Ｎ’、各粒子の質量ｍ’、減衰係数ｃ’が、
【数９】
0006472346.tif 000051
となるように、くりこみ変換を行うプログラムを記録したコンピュータ読み取り可能な記録媒体。
【請求項６】
各粒子に働く力が、粒子間の相互作用ポテンシャルφに依存するポテンシャル依存項と、エネルギの散逸に依存する散逸項とで表されるシミュレーション対象の粉粒体Ｓを定義する物理量を取得する手順と、
前記粉粒体Ｓの各粒子の運動エネルギと、前記相互作用ポテンシャルφに基づくポテンシャルエネルギとで表されるハミルトニアンの形が不変となるように、前記粉粒体Ｓの各粒子の運動方程式の前記ポテンシャル依存項に含まれる前記物理量をくりこみ変換する手順と、
前記物理量のくりこみ変換による前記ポテンシャル依存項の変化率が、前記散逸項の変化率と等しくなるように、前記散逸項に含まれる物理量をくりこみ変換する手順と、
くりこまれた粉粒体Ｓ’の各粒子の位置ベクトル及び運動量ベクトルに初期値を与えて、ある時間刻み幅で、くりこまれた前記物理量に基づいて、個々の粒子の運動方程式を数値積分することにより、くりこまれた粉粒体Ｓ’の時間発展を求める手順と
をコンピュータに実行させるプログラムを記録したコンピュータ読み取り可能な記録媒体であって、
前記相互作用ポテンシャルφが、粒子間距離に依存する無次元関数ｆ、エネルギの次元をもつ相互作用係数ε、前記粉粒体Ｓを特徴づけるパラメータｒ０、σ、及び粒子間距離ｒを用いて、
【数１０】
0006472346.tif 000052
と表すことができ、
前記散逸項が、減衰係数ｃと、２つの粒子の相対速度との積で表され、前記粉粒体Ｓの粒子数をＮ、各粒子の質量をｍで表し、前記粉粒体Ｓが含まれる空間の次元数をｄで表し、くりこみ回数に依存するくりこみ因子をαで表したとき、
繰り込まれた粉粒体Ｓ’の粒子数Ｎ’、各粒子の質量ｍ’、相互作用係数ε’、粉粒体Ｓ’を特徴づけるパラメータσ’、及び減衰係数ｃ’が、
【数１１】
0006472346.tif 000053
となるように、くりこみ変換を行い、くりこまれた粉粒体Ｓ’の粒子間の相互作用ポテンシャルφ’を、
【数１２】
0006472346.tif 000054
として、くりこまれた粉粒体Ｓ’の個々の粒子の運動方程式を数値積分するプログラムを記録したコンピュータ読み取り可能な記録媒体。
【請求項７】
各粒子に働く力が、粒子間の相互作用ポテンシャルφに依存するポテンシャル依存項と、エネルギの散逸に依存する散逸項とで表されるシミュレーション対象の粉粒体Ｓの時間発展をシミュレーションする装置であって、
入力装置と、
出力装置と、
処理ユニットと
を有し、
前記処理ユニットは、
前記入力装置を通して入力された前記粉粒体Ｓを定義する物理量を、前記粉粒体Ｓの各粒子の運動エネルギと、相互作用ポテンシャルφに基づくポテンシャルエネルギとで表されるハミルトニアンの形が不変となるように、前記粉粒体Ｓの各粒子の運動方程式のポテンシャル依存項に含まれる前記物理量をくりこみ変換し、
前記物理量のくりこみ変換による前記ポテンシャル依存項の変化率が、前記散逸項の変化率と等しくなるように、前記散逸項に含まれる物理量をくりこみ変換し、
くりこまれた前記物理量に基づいて、個々の粒子の運動方程式を、ある時間刻み幅で数値積分することにより、くりこまれた粉粒体Ｓ’の時間発展を求め、
求められた粉粒体Ｓ’の時間発展を、前記出力装置に出力するシミュレーション装置であって、
前記相互作用ポテンシャルφが、粒子間距離に依存する無次元関数ｆと、エネルギの次元をもつ相互作用係数εとの積εｆで表され、前記散逸項が、減衰係数ｃと、
２つの粒子の相対速度との積で表され、前記粉粒体Ｓの粒子数をＮ、各粒子の質量をｍで表し、前記粉粒体Ｓが含まれる空間の次元数をｄで表し、くりこみ回数に依存するくりこみ因子をαで表したとき、
繰り込まれた粉粒体Ｓ’の粒子数Ｎ’、各粒子の質量ｍ’、減衰係数ｃ’が、
【数１３】
0006472346.tif 000055
となるように、くりこみ変換を行うシミュレーション装置。
【請求項８】
各粒子に働く力が、粒子間の相互作用ポテンシャルφに依存するポテンシャル依存項と、エネルギの散逸に依存する散逸項とで表されるシミュレーション対象の粉粒体Ｓの時間発展をシミュレーションする装置であって、
入力装置と、
出力装置と、
処理ユニットと
を有し、
前記処理ユニットは、
前記入力装置を通して入力された前記粉粒体Ｓを定義する物理量を、前記粉粒体Ｓの各粒子の運動エネルギと、相互作用ポテンシャルφに基づくポテンシャルエネルギとで表されるハミルトニアンの形が不変となるように、前記粉粒体Ｓの各粒子の運動方程式のポテンシャル依存項に含まれる前記物理量をくりこみ変換し、
前記物理量のくりこみ変換による前記ポテンシャル依存項の変化率が、前記散逸項の変化率と等しくなるように、前記散逸項に含まれる物理量をくりこみ変換し、
くりこまれた前記物理量に基づいて、個々の粒子の運動方程式を、ある時間刻み幅で数値積分することにより、くりこまれた粉粒体Ｓ’の時間発展を求め、
求められた粉粒体Ｓ’の時間発展を、前記出力装置に出力するシミュレーション装置であって、
前記相互作用ポテンシャルφが、粒子間距離に依存する無次元関数ｆ、エネルギの次元をもつ相互作用係数ε、前記粉粒体Ｓを特徴づけるパラメータｒ０、σ、及び粒子間距離ｒを用いて、
【数１４】
0006472346.tif 000056
と表すことができ、
前記散逸項が、減衰係数ｃと、２つの粒子の相対速度との積で表され、前記粉粒体Ｓの粒子数をＮ、各粒子の質量をｍで表し、前記粉粒体Ｓが含まれる空間の次元数をｄで表し、くりこみ回数に依存するくりこみ因子をαで表したとき、
繰り込まれた粉粒体Ｓ’の粒子数Ｎ’、各粒子の質量ｍ’、相互作用係数ε’、粉粒体Ｓ’を特徴づけるパラメータσ’、及び減衰係数ｃ’が、
【数１５】
0006472346.tif 000057
となるように、くりこみ変換を行い、くりこまれた粉粒体Ｓ’の粒子間の相互作用ポテンシャルφ’を、
【数１６】
0006472346.tif 000058
として、くりこまれた粉粒体Ｓ’の個々の粒子の運動方程式を数値積分するシミュレーション装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、粉粒体の時間変化をシミュレーションする方法及び装置に関する。
【背景技術】
【０００２】
従来、粉粒体の時間変化のシミュレーションに、個別要素法が多く適用されていた。また、粉粒体の流動化について、粒子動力学を適用したシミュレーションが行われている。下記の非特許文献１に、薄い粉体層に振動を与えたときの粉体層の挙動を粒子動力学シミュレーションによって求める方法が示されている。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特許第４６６６３５７号公報
【非特許文献】
【０００４】
【非特許文献１】Physical Review Letters Vol. 77 (1996) pp.4166-4169, Keiko M. Aoki et al., "Spontaneous Wave Pattern Formation in Vibrated Granular Materials"
【発明の概要】
【発明が解決しようとする課題】
【０００５】
シミュレーション対象の粉粒体の粒子数が多くなると、巨視的な挙動をシミュレーションする際の計算時間が膨大になってしまう。上記特許文献１に、ハミルトニアンにくりこみ群変換を行うことで、ハミルトニアンの形を不変に保ち、粒子数を圧縮することができるくりこみ群分子動力学法が開示されている。くりこみ群変換に際し、ハミルトニアンの形が不変であるということは、繰り込まれた粉粒体の分解能は低下するが、繰り込まれた粉粒体の挙動は、元の粉粒体の挙動と相似になることを意味する。
【０００６】
上記特許文献１に示された分子動力学シミュレーションに適用される粒子間に働く力は、粒子間の相互作用ポテンシャルに依存する。これに対し、粒子動力学シミュレーションにおいては、粒子間に働く力は、粒子間相互作用ポテンシャルに依存する項と、エネルギの散逸を表す項とで表される。粒子間に働く力が、エネルギの散逸を表す項を含むため、粒子動力学シミュレーションに、上記特許文献１のくりこみ群変換を適用したのでは、粉粒体の挙動の相似性を保つことはできない。
【０００７】
本発明の目的は、粒子動力学シミュレーションの計算時間を短縮することが可能な粉粒体のシミュレーション方法を提供することである。本発明の他の目的は、粒子動力学シミュレーションの計算時間を短縮することが可能な粉粒体のシミュレーション装置を提供することである。
【課題を解決するための手段】
【０００８】
本発明の一観点によると、
各粒子に働く力が、粒子間の相互作用ポテンシャルφに依存するポテンシャル依存項と、エネルギの散逸に依存する散逸項とで表されるシミュレーション対象の粉粒体Ｓを定義し、
前記粉粒体Ｓの各粒子の運動エネルギと、前記相互作用ポテンシャルφに基づくポテンシャルエネルギとで表されるハミルトニアンの形が不変となるように、前記粉粒体Ｓの各粒子の運動方程式の前記ポテンシャル依存項に含まれる物理量をくりこみ変換し、
前記物理量のくりこみ変換による前記ポテンシャル依存項の変化率が、前記散逸項の変化率と等しくなるように、前記散逸項に含まれる物理量をくりこみ変換し、
くりこまれた粉粒体Ｓ’の各粒子の位置ベクトル及び運動量ベクトルに初期値を与えて、ある時間刻み幅で、くりこまれた前記物理量に基づいて、個々の粒子の運動方程式を数値積分することにより、くりこまれた粉粒体Ｓ’の時間発展を求める粉粒体のシミュレーション方法であって、
前記相互作用ポテンシャルφが、粒子間距離に依存する無次元関数ｆと、エネルギの次元をもつ相互作用係数εとの積εｆで表され、前記散逸項が、減衰係数ｃと、
２つの粒子の相対速度との積で表され、前記粉粒体Ｓの粒子数をＮ、各粒子の質量をｍで表し、前記粉粒体Ｓが含まれる空間の次元数をｄで表し、くりこみ回数に依存するくりこみ因子をαで表したとき、
繰り込まれた粉粒体Ｓ’の粒子数Ｎ’、各粒子の質量ｍ’、減衰係数ｃ’が、
【数３８】
0006472346.tif 000002
となるように、くりこみ変換を行う粉粒体のシミュレーション方法が提供される。
【０００９】
本発明の他の観点によると、
各粒子に働く力が、粒子間の相互作用ポテンシャルφに依存するポテンシャル依存項と、エネルギの散逸に依存する散逸項とで表されるシミュレーション対象の粉粒体Ｓの時間発展をシミュレーションする装置であって、
入力装置と、
出力装置と、
処理ユニットと
を有し、
前記処理ユニットは、
前記入力装置を通して入力された前記粉粒体Ｓを定義する物理量を、前記粉粒体Ｓの各粒子の運動エネルギと、相互作用ポテンシャルφに基づくポテンシャルエネルギとで表されるハミルトニアンの形が不変となるように、前記粉粒体Ｓの各粒子の運動方程式のポテンシャル依存項に含まれる前記物理量をくりこみ変換し、
前記物理量のくりこみ変換による前記ポテンシャル依存項の変化率が、前記散逸項の変化率と等しくなるように、前記散逸項に含まれる物理量をくりこみ変換し、
くりこまれた前記物理量に基づいて、個々の粒子の運動方程式を、ある時間刻み幅で数値積分することにより、くりこまれた粉粒体Ｓ’の時間発展を求め、
求められた粉粒体Ｓ’の時間発展を、前記出力装置に出力するシミュレーション装置であって、
前記相互作用ポテンシャルφが、粒子間距離に依存する無次元関数ｆ、エネルギの次元をもつ相互作用係数ε、前記粉粒体Ｓを特徴づけるパラメータｒ０、σ、及び粒子間距離ｒを用いて、
【数３９】
0006472346.tif 000003
と表すことができ、
前記散逸項が、減衰係数ｃと、２つの粒子の相対速度との積で表され、前記粉粒体Ｓの粒子数をＮ、各粒子の質量をｍで表し、前記粉粒体Ｓが含まれる空間の次元数をｄで表し、くりこみ回数に依存するくりこみ因子をαで表したとき、
繰り込まれた粉粒体Ｓ’の粒子数Ｎ’、各粒子の質量ｍ’、相互作用係数ε’、粉粒体Ｓ’を特徴づけるパラメータσ’、及び減衰係数ｃ’が、
【数４０】
0006472346.tif 000004
となるように、くりこみ変換を行い、くりこまれた粉粒体Ｓ’の粒子間の相互作用ポテンシャルφ’を、
【数４１】
0006472346.tif 000005
として、くりこまれた粉粒体Ｓ’の個々の粒子の運動方程式を数値積分するシミュレーション装置が提供される。
【発明の効果】
【００１０】
物理量のくりこみ変換によるポテンシャル依存項の変化率が、散逸項の変化率と等しくなるように、散逸項に含まれる物理量のくりこみ変換を行うため、散逸項を含む粉粒体の時間発展をシミュレーションすることができる。このとき、粉粒体Ｓの各粒子の運動エネルギと、粉粒体Ｓの相互作用ポテンシャルφに基づくポテンシャルエネルギとで表されるハミルトニアンの形が、くりこみ変換によって不変であるため、シミュレーション対象の粉粒体Ｓの時間発展と、くりこまれた粉粒体Ｓ’の時間発展とが、相似の関係を持つ。
【図面の簡単な説明】
【００１１】
【図１】図１は、実施例によるシミュレーション方法のフローチャートである。
【図２】図２は、実施例によるシミュレーション方法を実行するシミュレーション装置（コンピュータ）、及びこのシミュレーション装置にシミュレーション方法を実行させるプログラムを記録した記録媒体のブロック図である。
【図３】図３Ａは、シミュレーション対象の粉粒体Ｓの初期状態を示す模式図であり、図３Ｂは、くりこみ変換を行わない粉粒体Ｓの時間展開をシミュレーションした結果を示す模式図であり、図３Ｃは、くりこまれた粉粒体Ｓ’の時間展開をシミュレーションした結果を示す模式図である。
【図４】図４は、粒子ｉ、粒子ｊ、及び粒子ｋがこの順番に一次元状に並んで配置されている状態を示す模式図である。
【図５】図５は、相互作用ポテンシャルφがレナードジョーンズ型ポテンシャルである場合の式（２１）の算出結果、及び式（２２）の数値積分結果を示すグラフである。
【図６】図６は、相互作用ポテンシャルφがモース型ポテンシャルである場合の式（２１）の算出結果、及び式（２２）の数値積分結果を示すグラフである。
【図７】図７Ａ、図７Ｂ、及び図７Ｃは、二次元格子を一次元格子に帰着させるポテンシャルムービングの方法について説明する模式図である。
【発明を実施するための形態】
【００１２】
実施例によるシミュレーション方法を説明する前に、実施例によるシミュレーション方法で適用される粒子動力学について簡単に説明する。粒子数Ｎのシミュレーション対象の粉粒体Ｓの粒子ｊに働く力Ｆｊは、以下の式で表される。
【数１】
0006472346.tif 000006
ここで、φは、粒子間の相互作用ポテンシャルを表し、ベクトルｑｊは、粒子ｊの位置ベクトル（座標）を表し、ベクトルｖｊは、粒子ｊの速度ベクトルを表し、ｃは、減衰係数を表す。Σは、粒子ｊと接触している粒子について合計することを意味する。式（１）の右辺の大かっこ内の第１項は、粒子間の相互作用ポテンシャルφに依存するポテンシャル依存項であり、第２項は、エネルギの散逸に依存する散逸項である。散逸項は、減衰係数ｃと、相互作用する２つの粒子の相対速度との積で表される。
【００１３】
上述の式（１）から粒子ｊの運動方程式が得られる。個々の粒子に対して、運動方程式を数値積分することにより、個々の粒子の位置ベクトルｑｊ、及び運動量ベクトルｐｊの時間変化が求まる。
【００１４】
相互作用ポテンシャルφは、粒子間距離ｒに依存する無次元関数ｆと、エネルギの次元をもつ相互作用係数εとの積εｆで表される。相互作用ポテンシャルφは、一例として、以下の式で与えられる。
【数２】
0006472346.tif 000007
ここで、ｒは粒子間の距離を表し、σは粉粒体Ｓを特徴づけるパラメータである。式（２）のｒ＜ｒ０のときの相互作用ポテンシャルφの大かっこ内の第１項及び第２項は、レナードジョーンズ型ポテンシャルに相当する。第３項の定数１／４は、ｒ＝ｒ０において、相互作用ポテンシャルφを滑らかに接続するためのものである。
【００１５】
次に、図１〜図３Ｃを参照して、実施例によるシミュレーション方法について説明する。
図１に、本実施例によるシミュレーション方法のフローチャートを示す。ステップ１０において、シミュレーション対象の粉粒体Ｓを定義する。具体的には、粉粒体Ｓを特徴づける物理量である粒子数Ｎ、各粒子の質量ｍ、相互作用ポテンシャルφ、減衰係数ｃを決定する。
【００１６】
ステップ１１において、粉粒体Ｓの各粒子の運動エネルギと、相互作用ポテンシャルφに基づくポテンシャルエネルギとで表されるハミルトニアンの形が不変となるように、粉粒体Ｓの各粒子の運動方程式のポテンシャル依存項に含まれる物理量をくりこみ変換する。このくりこみ変換則については、特許第４６６６３５７号公報に説明されている。ハミルトニアンの形が不変であるということは、くりこみ前の粉粒体と、くりこまれた粉粒体との運動が、相互に相似であることを意味する。
【００１７】
一例として、くりこまれた粉粒体Ｓ’の粒子数Ｎ’、各粒子の質量ｍ’、各粒子の位置ベクトルｑ’、運動量ベクトルｐ’は、以下のくりこみ変換則によって求まる。
【数３】
0006472346.tif 000008
ここで、ｄは粉粒体Ｓが含まれる空間の次元数を表し、αはくりこみ回数に依存するくりこみ因子を表す。くりこみ回数をｎで表すと、くりこみ因子αは以下の式で表される。
【数４】
0006472346.tif 000009
【００１８】
式（３）のくりこみ変換則によってくりこまれた粉粒体Ｓ’に適用される相互作用ポテンシャルφは、シミュレーション対象の粉粒体Ｓの相互作用ポテンシャルφと同一である。例えば、くりこみ変換の前後において、式（２）に含まれる相互作用係数ε、パラメータσは不変である。相互作用ポテンシャルφが不変であるということは、元の粉粒体Ｓの安定状態における粒子間距離と、くりこまれた粉粒体Ｓ’の安定状態における粒子間距離とが等しいことを意味する。このため、くりこまれた粉粒体Ｓ’の巨視的な寸法が、元の粉粒体Ｓの巨視的な寸法の１／α倍に縮小される。
【００１９】
次に、ステップ１２において、物理量のくりこみ変換によるポテンシャル依存項の変化率と、散逸項の変化率とが等しくなるように、散逸項に含まれる物理量、具体的には減衰係数ｃ及び速度ｖのくりこみ変換を行う。以下、散逸項に含まれる物理量のくりこみ変換則について説明する。
【００２０】
くりこみ変換によって相互作用ポテンシャルφは不変であるため、運動方程式のポテンシャル依存項も不変である。すなわち、物理量のくりこみ変換によるポテンシャル依存項の変化率は１である。式（１）に示すように、運動方程式の散逸項は、減衰係数ｃと速度ｖとの積で表される。式（３）に示すように、くりこみ変換によって、運動量ベクトルｐの大きさがα倍になり、質量ｍがα２倍になる。速度ｖは、ｐ／ｍに等しいため、速度ｖは、くりこみ変換によって１／α倍になる。減衰係数ｃと速度ｖとの積で表される散逸項が、くりこみ変換によって不変になるためには、粉粒体Ｓの減衰係数ｃと、くりこまれた粉粒体Ｓ’の減衰係数ｃ’との間に、以下のくりこみ変換則を適用すればよい。
【数５】
0006472346.tif 000010
【００２１】
次に、ステップ１３において、ステップ１１及びステップ１２でくりこまれた物理量に基づいて、くりこまれた粉粒体Ｓ’の時間発展をシミュレーションする。具体的には、粉粒体Ｓ’の各粒子の位置ベクトルｑ及び運動量ベクトルｐの初期値を設定する。くりこまれた物理量に基づいて、ある時間刻み幅で個々の粒子の運動方程式を数値積分する。この数値積分により、個々の粒子の位置ベクトルｑ及び運動量ベクトルｐの時間変化が求まる。
【００２２】
次に、ステップ１４において、シミュレーション結果を出力する。シミュレーション結果は、例えば各時刻における粒子の位置を、画像として表示することができる。または、粉粒体Ｓ’の巨視的な外形を画像として表示することも可能である。
【００２３】
図２に、実施例によるシミュレーション方法を実行するシミュレーション装置（コンピュータ）、及びこのシミュレーション装置に、図１の各ステップを実行させるプログラムを記録した記録媒体のブロック図を示す。
【００２４】
シミュレーション装置２０は、読取装置２１、処理ユニット２２、入力装置２３、記憶装置２４、及び出力装置２５を含む。コンピュータ読み取り可能な記録媒体３０に、図１に示したシミュレーション方法をシミュレーション装置２０に実行させるためのプログラムが記録されている。
【００２５】
読取装置２１が、記録媒体３０に記録されているプログラムを読み取る。読み取られたプログラムは、例えば記憶装置２４に格納され、処理ユニット２２によって実行される。以下、図１に示した各ステップ１０〜１４における具体的な手順について説明する。
【００２６】
ステップ１０（図１）において、シミュレーション対象の粉粒体Ｓを定義する物理量が、入力装置２３を通して取得される。これにより、シミュレーション対象の粉粒体Ｓが定義される。
【００２７】
ステップ１１（図１）において、くりこみ変換の回数ｎが、入力装置２３を通して取得される。処理ユニット２２は、式（４）を用いてくりこみ因子αを算出する。算出されたくりこみ因子αを用いて、運動方程式のポテンシャル依存項に含まれる物理量のくりこみ変換を行う。このくりこみ変換は、式（３）に基づいて行うことができる。ステップ１２において、運動方程式の散逸項に含まれる物理量のくりこみ変換を行う。このくりこみ変換は、式（５）を用いて行うことができる。
【００２８】
ステップ１３において、処理ユニット２２が、くりこまれた粉粒体Ｓ’についてシミュレーションを実行する。ステップ１４において、処理ユニット２２が、シミュレーション結果を出力装置２５に出力する。
【００２９】
次に、図３Ａ〜図３Ｃを参照して、上記実施例による方法でシミュレーションを行った結果と、くりこみ変換を行わないでシミュレーションを行った結果とを比較して説明する。
【００３０】
図３Ａに、シミュレーション対象の粉粒体Ｓの初期状態を示す。仮想的な壁４０の中に、粉粒体Ｓが収納されている。仮想的な壁４０は、鉛直方向に延びる直方体の側面を構成する。粉粒体Ｓの各粒子の質量ｍを６．３８×１０２１／ＮＡ［ｇ］、相互作用係数εを６．８４７４×１０２２［Ｋ］、相互作用ポテンシャルφを決定するパラメータσを３．０×１０−３［ｍ］、減衰係数ｃを４．２４５３×１０４［ｋｇ／ｓ］とした。ここで、ＮＡはアボガドロ定数である。空間の次元数ｄを３とした。
【００３１】
シミュレーション対象の粉粒体Ｓの各粒子の運動方程式に含まれる物理量に対して１回のくりこみ変換を行った。すなわち、くりこみ回数ｎが１であり、くりこみ因子αが２である。シミュレーション対象の粉粒体Ｓの粒子数が約１７万個であり、くりこまれた粉粒体Ｓ’の粒子数が約２．１万個である。シミュレーションの開始時刻において、仮想的な壁４０を取り外し、その後の粉粒体Ｓ及びＳ’の時間展開をシミュレーションにより求めた。粉粒体Ｓ及びＳ’が定常状態になった時点でシミュレーションを終了した。
【００３２】
図３Ｂに、くりこみ変換を行わない粉粒体Ｓの時間展開をシミュレーションした結果を示し、図３Ｃに、くりこまれた粉粒体Ｓ’の時間展開をシミュレーションした結果を示す。くりこまれた粉粒体Ｓ’の巨視的な寸法は、元の粉粒体Ｓの巨視的な寸法の１／２に縮小され、体積は１／８に縮小されている。図３Ｂと図３Ｃとを比較すると、両者の巨視的な形状はほぼ相似であることがわかる。従って、くりこまれた粉粒体Ｓ’の時間展開のシミュレーション結果から、元の粉粒体Ｓの時間展開を推測することができる。
【００３３】
図３Ｂに示した粉粒体Ｓのシミュレーションに要した計算時間は約１２１時間であった。これに対し、図３Ｃに示したくりこまれた粉粒体Ｓ’のシミュレーションに要した計算時間は約２６時間であった。このように、シミュレーション対象の粉粒体Ｓに対してくりこみ変換を行うことにより、シミュレーションに要する計算時間を短縮することができる。くりこみ回数を増やすと、計算時間をより短縮することが可能である。
【００３４】
次に、図４〜図７Ｃを参照して、他の実施例によるシミュレーション方法について説明する。以下、図１〜図３Ｃに示した実施例との相違点について説明し、共通の構成については説明を省略する。本実施例においても、粉粒体Ｓの粒子ｊに加わる力Ｆｊは、上述の式（１）で表される。
【００３５】
図１〜図３Ｃに示した実施例では、元の粉粒体Ｓとくりこまれた粉粒体Ｓ’とで、相互作用ポテンシャルφが同一であり、くりこまれた粉粒体Ｓ’の巨視的な寸法が元の粉粒体Ｓの巨視的な寸法の１／α倍に縮小された。以下に説明する実施例では、くりこまれた粉粒体Ｓ’の巨視的な寸法と元の粉粒体Ｓの巨視的な寸法とが同一であり、相互作用ポテンシャルφが異なる。ただし、本実施例においても、粉粒体Ｓの各粒子の運動エネルギと、相互作用ポテンシャルφに基づくポテンシャルエネルギとで表されるハミルトニアンの形は、くりこみ変換によって不変である。
【００３６】
以下、本実施例で適用されるくりこみ変換則について説明する。シミュレーション対象の粉粒体Ｓに適用される相互作用ポテンシャルφが、粒子間距離ｒに依存する無次元関数ｆと、相互作用係数εとを用いて以下の式で定義される。
【数６】
0006472346.tif 000011
【００３７】
上述の式（２）は、下記のように書き直すことができる。
【数７】
0006472346.tif 000012
従って、式（２）で表される相互作用ポテンシャルφも、式（６）で示す相互作用ポテンシャルφに含まれる。このため、相互作用ポテンシャルが式（２）で表される粉粒体Ｓも、本実施例で取り扱うことができる。
【００３８】
ステップ１０（図１）において定義される粉粒体Ｓを特徴づける物理量として、図１〜図３Ｃに示した実施例で粉粒体Ｓが定義された物理量に加えて、式（６）のパラメータｒ０が定義される。
【００３９】
ステップ１１（図１）において、下記のくりこみ変換則により、粉粒体Ｓを定義する物理量、である粒子数Ｎ、粒子の質量ｍ、相互作用係数ε、無次元関数ｆを規定するパラメータｒ０及びσが、それぞれ粉粒体Ｓ’を定義する物理量である粒子数Ｎ’、粒子の質量ｍ’、相互作用係数ε’、無次元関数ｆを規定するパラメータｒ０’及びσ’に変換される。
【数８】
0006472346.tif 000013
【００４０】
くりこまれた粉粒体Ｓ’の相互作用ポテンシャルφ’は、以下の式で表すことができる。無次元関数ｆの形は不変である。
【数９】
0006472346.tif 000014
【００４１】
くりこまれた粉粒体Ｓ’のハミルトニアンＨ’は、後に詳細に説明するように、以下の式で表すことができる。ハミルトニアンの形は、くりこみ変換によって不変である。
【数１０】
0006472346.tif 000015
【００４２】
式（８）のくりこみ変換測を適用することにより、粒子の個数が１／αｄ倍になり、粒子の質量がαｄ倍になる。このため、粉粒体Ｓの全体の質量と、粉粒体Ｓ’の全体の質量とは同一である。さらに、パラメータｒ０、すなわち安定状態における粒子間距離がα倍になる。このため、粉粒体Ｓとくりこまれた粉粒体Ｓ’との巨視的な寸法も同一である。くりこみ変換処理の前後で、粉粒体の全質量及び寸法が不変であるため、粉粒体の密度も不変である。
【００４３】
次に、ステップ１２（図１）で行われる散逸項に含まれる物理量のくりこみ変換について説明する。運動方程式のポテンシャル依存項は、以下のように変形される。
【数１１】
0006472346.tif 000016
ここで、Ｄｆは以下の式で定義される。
【数１２】
0006472346.tif 000017
【００４４】
くりこまれた粉粒体Ｓ’に含まれる各粒子の運動方程式のポテンシャル依存項は、以下のように変形される。
【数１３】
0006472346.tif 000018
【００４５】
式（１１）と式（１３）とから、以下の式が導出される。
【数１４】
0006472346.tif 000019
式（１４）から、くりこみ変換によって、ポテンシャル依存項がα（ｄ−１）倍になることがわかる。
【００４６】
式（８）の質量ｍのくりこみ変換則、及び式（１０）の運動量のくりこみ変換則から、くりこみ変換によって速度ｖは不変であることがわかる。式（１）の散逸項に含まれる速度ｖが不変であるから、散逸項をα（ｄ−１）倍にするためには、粉粒体Ｓの減衰係数ｃから、粉粒体Ｓ’の減衰係数ｃ’へのくりこみ変換則を以下の式で定義すればよい。
【数１５】
0006472346.tif 000020
【００４７】
本実施例によるステップ１３及びステップ１４（図１）の処理は、図１〜図３Ｃに示した処理と共通である。
【００４８】
次に、図４〜図７を参照して、くりこまれた粉粒体Ｓ’のハミルトニアンＨ’（式（１０））の導出について説明する。粉粒体Ｓに対する分配関数Ｚ（β）の積分の一部を実行し、ハミルトニアンを粗視化することにより、くりこまれた粉粒体Ｓ’のハミルトニアンＨ’が得られる。
【００４９】
粒子数が一定の正準集団（Canonical ensemble）に対する分配関数Ｚ（β）は、以下の式で表される。
【数１６】
0006472346.tif 000021
ここで、ｄΓＮは位相空間内の体積要素であり、下記の式で表される。
【数１７】
0006472346.tif 000022
ここで、ｈはプランク定数である。ＷＮは、すべての状態の量子論的な和と、位相空間に亘る積分とが一致するように決められる。
【００５０】
まず、粒子間の相互作用ポテンシャルの粗視化について説明し、次に運動エネルギの粗視化について説明する。続いて、相互作用ポテンシャルの粗視化、及び運動エネルギの粗視化を踏まえて、くりこみ変換測を定義する。
【００５１】
［粒子間の相互作用ポテンシャルの粗視化］
まず、一次元鎖状に配置された粉粒体における相互作用ポテンシャルの粗視化について説明する。その後、単純立方格子状に配置された粉粒体における相互作用ポテンシャルについて説明する。
【００５２】
図４に示すように、粒子ｉ、粒子ｊ、及び粒子ｋがこの順番に一次元状に並んで配置されている。粒子ｊが関与する相互作用を書き出し、粒子ｉと粒子ｋとの中間にある粒子ｊの位置座標について積分を実行することにより、相互作用ポテンシャルの粗視化を行うことができる。まず、第二近接以上の粒子からの寄与を取り込むために、ポテンシャルムービング（Potential Moving）の方法を用いることができる。ポテンシャルムービングの方法については、Leo P. Kadanoff, STATISTICAL PHYSICS Static, Dynamics and Renormalization, Chap. 14, World Scientific (1999)に解説されている。
【００５３】
第二近接以上の粒子からの寄与を取り込んだ相互作用ポテンシャルφティルダは、下記の式で表すことができる。
【数１８】
0006472346.tif 000023
ここで、ａは平衡状態における粒子間距離を表す。平衡状態における粒子間距離ａは、相互作用ポテンシャルφが最小となる距離であるｒ０に等しいと近似することができる。
【００５４】
複数の粒子が一次元状に配置されているため、粒子ｊの位置ベクトルｑｊは、一次元座標で表すことができる。粒子ｊの位置をｑｊで表すと、粒子ｊに対して、最近接の粒子によって作られるケージポテンシャルは、下記の式で表される。
【数１９】
0006472346.tif 000024
【００５５】
位置ベクトルｑｊについて積分を実行すれば、式（１９）を用いて下記の式が得られる。
【数２０】
0006472346.tif 000025
ここで、ｒａは粒子の直径を表し、ｚ（ｑｉ−ｑｋ）及びＰ（ｑｉ−ｑｋ）は、下記の式で表される。
【数２１】
0006472346.tif 000026
【数２２】
0006472346.tif 000027
積分領域はケージポテンシャルの内部領域に制限される。
【００５６】
次に、ｚ（ｑｉ−ｑｋ）を具体的に計算する。相互作用ポテンシャルφがレナードジョーンズ型ポテンシャルの場合、φ（２ｎ）は、下記の式で表される。
【数２３】
0006472346.tif 000028
相互作用ポテンシャルφがモース型ポテンシャルの場合、φ（２ｎ）は、下記の式で表される。
【数２４】
0006472346.tif 000029
【００５７】
式（２３）または式（２４）を式（２２）に代入して数値積分を行う。式（２２）に式（２３）または式（２４）を代入する際には、式（１８）を用いる。数値積分において、式（２２）の積分範囲に現れているａは、近似的にｒ０に等しいと仮定する。
【００５８】
図５に、相互作用ポテンシャルφがレナードジョーンズ型ポテンシャルである場合の式（２１）の算出結果、及び式（２２）の数値積分結果を示す。図６に、相互作用ポテンシャルφがモース型ポテンシャルである場合の式（２１）の算出結果、及び式（２２）の数値積分結果を示す。粒子ｉ、粒子ｊ、粒子ｋ、粒子ｌ、及び粒子ｍがこの順番に一次元状に配置されている場合の粒子ｋの位置座標がｑｋで表される。図５及び図６の横軸は、ｑｋ／２を表し、縦軸は、Ｐ（ｑｉ−ｑｋ）Ｐ（ｑｋ−ｑｍ）、及びｚ（ｑｉ−ｑｋ）ｚ（ｑｋ−ｑｍ）を、対数目盛で表す。図５に至る数値計算において、ε／ｋＢＴ＝２．０、ｒ０／σ＝１．１２とした。図６に至る数値計算において、ε／ｋＢＴ＝２．０、ｒ０／σ＝２．２４とした。
【００５９】
相互作用ポテンシャルφがレナードジョーンズ型ポテンシャルである場合、及びモース型ポテンシャルである場合のいずれにおいても、Ｐ（ｑｉ−ｑｋ）Ｐ（ｑｋ−ｑｍ）の変化に比べて、ｚ（ｑｉ−ｑｋ）ｚ（ｑｋ−ｑｍ）の変化が緩やかであることがわかる。このため、Ｐ（ｑｉ−ｑｋ）Ｐ（ｑｋ−ｑｍ）に対してｚ（ｑｉ−ｑｋ）ｚ（ｑｋ−ｑｍ）は、ほぼ定数であると近似できる。
【００６０】
粒子ｋが位置座標ｑｋに存在する確率ｐ（ｑｋ）は、以下のように近似することができる。
【数２５】
0006472346.tif 000030
【００６１】
従って、下記の式が導出される。
【数２６】
0006472346.tif 000031
【００６２】
ここまでは、一次元状に複数の粒子が配置された粉粒体の相互作用ポテンシャルの粗視化について説明した。多次元の粉粒体の相互作用ポテンシャルは、ポテンシャルムービングの方法により実現できる。
【００６３】
図７Ａ〜図７Ｃを参照して、二次元格子を一次元格子に帰着させるポテンシャルムービングの方法について説明する。
【００６４】
図７Ａに示すように、二次元正方格子の格子点の位置に粒子が配置されている。最近接する粒子間の相互作用（ニアレストネイバーカップリング）を実線で示す。粒子が配列する１つの方向をｘ方向とし、それに直交する方向をｙ方向とする。
【００６５】
図７Ｂに示すように、ｘ方向に並んだ粒子のうち、１つ置きに配置された粒子（中空の円で表された粒子）の変位について積分を実行することを考える。積分の対象となる粒子（消去したい粒子）を被積分粒子と呼ぶこととする。
【００６６】
図７Ｃに示すように、ポテンシャルムービングの方法では、被積分粒子同士のニアレストネイバーカップリングを、ｘ方向に関して両隣に配置された粒子に振り分ける。振り分けられた相互作用が加算された粒子間相互作用（ダブルカップリング）を二重線で表す。二重線で表されたダブルカップリングは、元のニアレストネイバーカップリングの２倍の強さになる。この手法を用いて、二次元格子を一次元鎖に変換することができる。一次元鎖の粉粒体においては、図４を参照して説明した手法により相互作用ポテンシャルの粗視化を行うことができる。シミュレーション対象の粉粒体が三次元格子を構成する場合には、二次元格子を一次元鎖に変換する手順をｘ方向、ｙ方向、及びｚ方向の３方向に関して繰り返せばよい。このようにして、多次元格子を構成する粉粒体の粗視化を行うことができる。
【００６７】
次元数ｄの多次元格子を構成する粉粒体の相互作用ポテンシャルの粗視化は、下記の式で表される。
【数２７】
0006472346.tif 000032
ここで、＜ｉ，ｊ＞は、最近接格子間で和を取ることを意味する。
【００６８】
すべての相互作用の和に変更すると、下記の式が得られる。
【数２８】
0006472346.tif 000033
【００６９】
［運動エネルギの粗視化］
次に、運動エネルギの粗視化について説明する。運動エネルギについては、容易に積分を実行することができ、下記の式が導出される。
【数２９】
0006472346.tif 000034
式（２９）の導出に際し、下記の式が利用される。ここで、運動量ベクトルｐｊ２は、ベクトルの内積を意味する。
【数３０】
0006472346.tif 000035
【００７０】
［くりこみ変換測の導出］
次に、上述の相互作用ポテンシャルの粗視化、及び運動エネルギの粗視化から導出されるくりこみ変換測について説明する。
【００７１】
式（１６）に、式（２８）及び式（２９）を代入することにより、結果に影響しない係数を除いて下記の式が得られる。
【数３１】
0006472346.tif 000036
【００７２】
式（３１）から、粗視化されたハミルトニアン（くりこまれた粉粒体Ｓ’のハミルトニアン）Ｈ’は、以下の式で表される。
【数３２】
0006472346.tif 000037
【００７３】
ハミルトニアンを粗視化する際の結合定数のリストをＫで表す。結合定数のリストＫは下記のように表される。
【数３３】
0006472346.tif 000038
【００７４】
くりこみ変換Ｒを以下のように定義する。
【数３４】
0006472346.tif 000039
【００７５】
ｎ回のくりこみ変換実行後の結合係数のリストＫｎは、下記の式で表される。
【数３５】
0006472346.tif 000040
【００７６】
従って、ｎ回のくりこみ変換が行われたハミルトニアンＨｎは、下記の式で表される。
【数３６】
0006472346.tif 000041
ここで、くりこまれた粉粒体Ｓ’における運動量ベクトルｐｊ’は、下記の式で表される。
【数３７】
0006472346.tif 000042
【００７７】
式（３６）から、式（８）に示したくりこみ変換測、及び式（１０）に示したくりこまれた粉粒体Ｓ’のハミルトニアンＨ’が導き出される。
【００７８】
本実施例においても、くりこまれた粉粒体Ｓ’の各粒子の運動エネルギと、相互作用ポテンシャルφに基づくポテンシャルエネルギとで表されるハミルトニアンの形が、くりこみ前の粉粒体Ｓのハミルトニアンの形と同一である。このため、本実施例においても、図１〜図３Ｃに示した実施例と同様に、くりこみ前の粉粒体Ｓと、くりこまれた粉粒体Ｓ’との運動が、相互に相似である。
【００７９】
以上実施例に沿って本発明を説明したが、本発明はこれらに制限されるものではない。例えば、種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。
【符号の説明】
【００８０】
１０、１１、１２、１３、１４  ステップ
２０  シミュレーション装置（コンピュータ）
２１  読取装置
２２  処理ユニット
２３  入力装置
２４  記憶装置
２５  出力装置
３０  記録媒体
４０  仮想的な壁
Ｈ’、Ｈｎ  ハミルトニアン
Ｎ、Ｎ’  粒子数
Ｓ  シミュレーション対象の粉粒体
Ｓ’  繰り込まれた粉粒体
ｃ、ｃ’  減衰係数
ｄ  次元数
ｆ  無次元関数
ｍ、ｍ’  質量
ｐ、ｐ’  運動量ベクトル
ｑ、ｑ’  位置ベクトル
ｒ  粒子間距離
ｖ  速度
α  くりこみ因子
ε、ε’  相互作用係数
φ、φ’  相互作用ポテンシャル
【図１】
0006472346.tif 000059
【図２】
0006472346.tif 000060
【図３】
0006472346.tif 000061
【図４】
0006472346.tif 000062
【図５】
0006472346.tif 000063
【図６】
0006472346.tif 000064
【図７】
0006472346.tif 000065