JP2016218767A - シミュレーション方法、シミュレーションプログラム、及びシミュレーション装置 - Google Patents

Abstract

【課題】粒子の速度分布が維持されるくりこみ変換を適用してシミュレーションを行う方法を提供する。【解決手段】くりこみ回数に依存するくりこみ因子αに基づいて、シミュレーション対象の粒子系Ｓに対してくりこみ変換処理を行う。くりこまれた粒子系Ｓ’について分子動力学計算を実行することにより、粒子系Ｓ’の粒子の位置ベクトル及び運動量ベクトルを求める。粒子系Ｓの粒子間の相互作用ポテンシャルφが、エネルギの次元を持つ相互作用係数ε、無次元関数ｆ、粒子を特徴づけるパラメータｒ０、σ、粒子間距離ｒを用いて、φ（ｒ）＝εｆ（（ｒ−ｒ０）／σ）と表される。粒子系Ｓの空間の次元をｄで表したとき、変換側Ｎ’＝Ｎ／αｄ、ｍ’＝ｍαｄ、ε’＝εαｄ、ｒｏ’＝αｒ０、σ’＝ασを適用し、粒子系Ｓ’の相互作用ポテンシャルφ’（ｒ）＝ε’ｆ（（ｒ−ｒ０’）／σ’）に基づいて分子動力学計算を実行する。【選択図】図１

Description

本発明は、くりこみ群を応用した分子動力学によるシミュレーション方法、シミュレーションプログラム、及びシミュレーション装置に関する。

分子動力学を用いたコンピュータシミュレーションが行われている。分子動力学では、シミュレーションの対象となる系を構成する粒子の運動方程式が数値的に解かれる。シミュレーション対象の系に含まれる粒子数が増加すれば、必要な計算量が増大する。現行のコンピュータの演算応力でシミュレーション可能な系の粒子数は、通常、数十万個程度である。

下記の特許文献１に、必要な計算量を減らすためにくりこみ変換の手法を用いたシミュレーション方法が開示されている。以下、特許文献１に開示されたくりこみ変換の手法について説明する。

シミュレーション対象の粒子系Ｓの粒子数をＮ、各粒子の質量をｍ、粒子間の相互作用ポテンシャルをφ（ｒ）で表す。ここで、ｒは、粒子間の距離である。相互作用ポテンシャルφ（ｒ）が、相互作用係数εと関数ｆ（ｒ）との積で表される。相互作用係数εは、相互作用の強さを表し、エネルギの次元を持つ。関数ｆ（ｒ）は、粒子間距離に対する依存性を表し、無次元化されている。

第１のくりこみ因子α、第２のくりこみ因子γ、及び第３のくりこみ因子δが決定される。第１のくりこみ因子αは１より大きい。第２のくりこみ因子γは、０以上であり、かつ空間の次元数ｄ以下である。第３のくりこみ因子δは、０以上である。くりこみの回数をｎで表したとき、第１のくりこみ因子αは、
α＝２^ｎ
と表される。

くりこみ手法を用いてくりこみ変換された粒子系Ｓ’の粒子数をＮ’、各粒子の質量をｍ’、相互作用係数をε’で表す。以下の変換式を用いて、くりこみ変換された粒子系Ｓ’の粒子数Ｎ’、質量ｍ’、及び相互作用係数ε’を求める。
Ｎ’＝Ｎ／α^ｄ
ｍ’＝ｍα^δ−γ
ε’＝εα^γ

くりこみ変換された粒子系Ｓ’について、分子動力学計算を行なう。分子動力学計算で求められた各粒子の位置ベクトルをｑ’、運動量ベクトルをｐ’で表す。元の粒子系Ｓの各粒子の位置ベクトルｑ、運動量ベクトルｐは、以下の式を用いて算出することができる。
ｑ＝ｑ’α
ｐ＝ｐ’／α^δ／２

特許第５２４１４６８号公報

くりこみ変換された粒子系Ｓ’のマクスウェル速度分布則ｆ_ｍａｘ（ｖ’）は、以下の式で表される。
ｆ_ｍａｘ（ｖ’）＝ｅｘｐ［（−ｍ’／２ｋ_ＢＴ）ｖ’^２］
ここで、ｖ’は、くりこみ変換された粒子系Ｓ’の各粒子の速度であり、ｋ_Ｂはボルツマン定数であり、Ｔは粒子系Ｓ’の温度である。

上述のマクスウェル速度分布則ｆ_ｍａｘ（ｖ’）は、元の粒子系Ｓの各粒子の質量ｍを用いて、以下の式に書き直される。
ｆ_ｍａｘ（ｖ’）＝ｅｘｐ［（−ｍ／２ｋ_ＢＴ）（ｖ’α^{（δ−γ）／２}）^２］
上述の式は、くりこみ変換された粒子系Ｓ’の粒子の速度分布の標準偏差が、元の粒子系Ｓの粒子の速度分布の標準偏差の１／α^{（δ−γ）／２}倍になっていることを意味する。すなわち、くりこみ変換された粒子系Ｓ’では、元の粒子系Ｓの粒子の速度分布が再現されていない。くりこみ回数ｎが大きくなると、くりこみ変換された粒子系Ｓ’の粒子の速度分布が、元の粒子系Ｓの粒子の速度分布から大きく乖離してしまう。

くりこみ変換された粒子系Ｓ’と元の粒子系Ｓとで速度分布が大きく異なると、蒸発量、液滴の生成、消滅等が正しく再現されなくなる。

本発明の目的は、粒子の速度分布が維持されるくりこみ変換を適用してシミュレーションを行う方法を提供することである。本発明の他の目的は、粒子の速度分布が維持されるくりこみ変換を適用してシミュレーションを行うコンピュータプログラムを提供することである。本発明のさらに他の目的は、粒子の速度分布が維持されるくりこみ変換を適用してシミュレーションを行う装置を提供することである。

本発明の一観点によると、
くりこみ回数に依存するくりこみ因子αに基づいて、複数の粒子からなるシミュレーション対象の粒子系Ｓに対してくりこみ変換処理を行う工程と、
くりこまれた粒子系Ｓ’について分子動力学計算を実行することにより、くりこまれた前記粒子系Ｓ’の粒子の位置ベクトル及び運動量ベクトルを求める工程と
を有し、
前記粒子系Ｓの粒子間の相互作用ポテンシャルφが、エネルギの次元を持つ相互作用係数ε、無次元関数ｆ、粒子を特徴づけるパラメータｒ_０、σ、粒子間距離ｒを用いて、

と表すことができ、
前記粒子系Ｓの空間の次元をｄで表したとき、下記の変換側

を適用し、くりこまれた前記粒子系Ｓ’の相互作用ポテンシャル

に基づいて前記分子動力学計算を実行するシミュレーション方法が提供される。

本発明の他の観点によると、上記シミュレーション方法を実行するコンピュータプログラムが提供される。本発明のさらに他の観点によると、上記コンピュータプログラムが記録された記録媒体が提供される。本発明のさらに他の観点によると、上記シミュレーション方法を実行するシミュレーション装置が提供される。

上述のシミュレーション方法を適用すると、くりまれた粒子系の速度分布が、元の粒子系の速度分布と同一になる。

図１は、実施例によるシミュレーション方法のフローチャートである。 図２は、一次元状に配置された粒子を示す概念図である。 図３は、相互作用ポテンシャルφがレナードジョーンズ型ポテンシャルである場合の式（２１）の算出結果、及び式（２２）の数値積分結果を示すグラフである。 図４は、相互作用ポテンシャルφがモース型ポテンシャルである場合の式（２１）の算出結果、及び式（２２）の数値積分結果を示すグラフである。 図５Ａ〜図５Ｃは、二次元格子状に配置された粒子、及び粒子間相互作用を示す概念図である。 図６Ａ〜図６Ｄは、実施例によるシミュレーション方法を用いたシミュレーション結果を示す図である。 図７Ａ〜図７Ｄは、比較例によるシミュレーション方法を用いたシミュレーション結果を示す図である。 図８は、実施例によるシミュレーション装置のブロック図である。

本願発明の実施例で適用される分子動力学（Molecular Dynamics）について簡単に説明する。Ｎ個の粒子（例えば原子）からなり、ハミルトニアンＨが下記の式で表される粒子系について考える。

ここで、
ｍは粒子の質量を表し、
φは粒子間の相互作用ポテンシャルを表し、
ベクトルｐ_ｊは粒子の運動量ベクトルを表し、
ベクトルｑ_ｊは粒子の位置ベクトル（位置座標）を表す。

ハミルトニアンＨを、ハミルトニアンの正準方程式に代入することにより、粒子ｊに対する下記の運動方程式が得られる。

分子動力学では、粒子系を構成する各粒子に対し、式（５）及び式（６）で表される運動方程式を数値的に積分して解くことにより、各時刻における各粒子の運動量ベクトルｐ_ｊ及び位置ベクトルｑ_ｊが求められる。数値積分には、Verlet法が用いられることが多い。Verlet法については、例えば、J.M.Thijssen,「Computational Physics」,Cambridge University Press (1999)の１７５頁に説明されている。分子動力学計算で得られた各粒子の運動量ベクトル及び位置ベクトルに基づいて、粒子系の種々の物理量を算出することができる。

次に、くりこみ郡の手法を応用した分子動力学（以下、くりこみ群分子動力学という。）について概念的に説明する。

くりこみ群分子動力学では、シミュレーション対象の粒子系Ｓを、粒子系Ｓの粒子数よりも少ない数の粒子で構成された粒子系Ｓ’（以下、くりこまれた粒子系Ｓ’という。）に対応づける。次に、くりこまれた粒子系Ｓ’に対して分子動力学計算を実行する。くりこまれた粒子系Ｓ’に対する計算結果を、シミュレーション対象の粒子系Ｓに対応付ける。これにより、シミュレーション対象の粒子系Ｓに対して分子動力学計算を直接的に実行する場合よりも計算量を少なくすることができる。シミュレーション対象の粒子系Ｓにおける物理量（例えば、粒子数、粒子の質量等）と、くりこまれた粒子系Ｓ’における物理量とを相互に対応付けるための変換則を、くりこみ変換測という。

図１に、実施例によるシミュレーション方法のフローチャートを示す。ステップＳ１において、シミュレーション対象の粒子系Ｓに対して、くりこみ変換測に基づいてくりこみ変換処理を実行することにより、くりこまれた粒子系Ｓ’を定義する。シミュレーション対象の粒子系Ｓにおいては、粒子間の相互作用ポテンシャルが下記の式で表される。

ここで、
ｒは粒子間の距離を表し、
ｆは無次元化された関数を表し、
ε、ｒ_０、σは、粒子（例えば、原子や分子）を特徴づけるパラメータである。εはエネルギの次元を持ち、相互作用係数と呼ばれる。ｒ_０は、相互作用ポテンシャルφが最小となる位置に相当する。平衡状態において、粒子間距離はｒ_０にほぼ等しい。

シミュレーション対象の粒子系Ｓの粒子が不活性原子の場合には、相互作用ポテンシャルφとして、レナードジョーンズ型ポテンシャルを適用することができる。レナードジョーンズ型ポテンシャルは、例えば下記の式で定義される。

式（８）は、下記の式のように変形することができる。

式（９）からわかるように、レナードジョーンズ型ポテンシャルは（ｒ−ｒ_０）／σの関数であり、式（７）の形に表すことができる。

シミュレーション対象の粒子系Ｓの粒子が金属原子の場合には、相互作用ポテンシャルφとして、モース型ポテンシャルを適用することができる。モース型ポテンシャルは、例えば下記の式で定義される。

くりこみ変換処理により、シミュレーション対象の粒子系Ｓの物理量Ｎ、ｍ、ε、ｒ_０、及びσが、それぞれ繰り込まれた粒子系Ｓ’の物理量Ｎ’、ｍ’、ε’、ｒ_０’、及びσ’に変換される。ステップＳ１のくりこみ変換処理では、下記のくりこみ変換測が適用される。

ここで、ｄはシミュレーション対象の粒子系Ｓが配置された空間の次元数を表す。αは繰り込み回数に依存するくりこみ因子である。くりこみ回数がｎ回のとき、くりこみ因子αは以下の式で表される。

くりこまれた粒子系Ｓ’の相互作用ポテンシャルφ’は、下記の式で表すことができる。

くりこまれた粒子系Ｓ’のハミルトニアンＨ’は、後に詳細に説明するように、以下の式で表すことができる。

上述のくりこみ変換測を適用することにより、粒子の個数が１／α^ｄ倍になり、粒子の質量がα^ｄ倍になる。このため、粒子系Ｓの全体の質量と、粒子系Ｓ’の全体の質量とは同一である。さらに、粒子間距離ｒ_０がα倍になる。このため、粒子系Ｓとくりこまれた粒子系Ｓ’との寸法は同一である。くりこみ変換処理の前後で、粒子系の全質量及び寸法が不変であるため、粒子系の密度も不変である。

次に、ステップＳ２において、シミュレーションの初期条件を設定する。初期条件には、各粒子の位置ベクトルｑ_ｊ及び運動量ベクトルｐ_ｊの初期値が含まれる。運動量ベクトルｐ_ｊは、くりこまれた粒子系Ｓ’の温度Ｔ’に基づいて設定される。シミュレーション対象の粒子系Ｓの温度をＴで表したとき、温度に関して下記のくりこみ変換測が適用される。

次に、ステップＳ３において、くりこまれた粒子系Ｓ’に関して分子動力学計算を実行する。具体的には、式（１３）で表されるくりこまれた粒子系Ｓ’のハミルトニアンＨ’を正準方程式に代入して運動方程式を得る。運動方程式は、下記の式で表される。

この運動方程式を、数値的に積分して解く。これにより、くりこまれた粒子系Ｓ’の各粒子の位置ベクトルｑ’及び運動量ベクトルｐ’の時刻歴が求まる。

くりこまれた粒子系Ｓ’の各粒子の位置ベクトルｑ’、運動量ベクトルｐ’と、シミュレーション対象の粒子系Ｓの位置ベクトルｑ、運動量ベクトルｐとは、下記の関係を有する。

ステップＳ４において、シミュレーション結果を出力する。例えば、位置ベクトルｑ’及び運動量ベクトルｐ’を、そのまま数値で出力してもよいし、位置ベクトルｑ’に基づいて、粒子系Ｓ’の複数の粒子の空間内の分布を画像として表示してもよい。その他に、シミュレーション結果に基づいてアクチュエータを駆動することにより、対象物に物理的作用を加える（オブザーバとする）ことも可能である。

次に、くりこまれた粒子系Ｓ’のハミルトニアンＨ’の導出について説明する。くりこまれた粒子系Ｓ’のハミルトニアンＨ’を得るためには、粒子系Ｓに対する分配関数Ｚ（β）の積分の一部を実行し、ハミルトニアンを粗視化することによりハミルトニアンＨ’が得られる。

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

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

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

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

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

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

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

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

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

積分変数ｑ_ｊについて積分を実行すれば、式（２１）を用いて下記の式が得られる。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

次に、上記実施例によるシミュレーション方法の優れた効果について説明する。くりこまれた粒子系Ｓ’におけるマクスウェルの速度分布則は、下記の式で表される。

式（４０）に、式（１１）のｍ’に関する式、及び式（１５）を代入することにより、下記の式が得られる。

式（４１）からわかるように、くりこまれた粒子系Ｓ’の速度分布は、くりこみ前の元の粒子系Ｓの速度分布と同一である。このため、粒子系Ｓの界面における粒子の挙動に関する現象、例えば蒸発量、液滴の生成と消滅等の再現性を高めることができる。

さらに、上記実施例によるシミュレーション方法では、粒子数Ｎが１／α^ｄ倍に変換され、平衡状態における粒子間距離ｒ_０がα倍に変換される。このため、くりこみ変換の前後において粒子系の寸法が不変である。また、粒子数Ｎが１／α^ｄ倍に変換され、粒子の質量がα^ｄ倍に変換される。このため、くりこみ変換の前後において粒子系の密度が不変である。

粒子系の寸法及び密度が不変であるため、上記実施例によるシミュレーション方法と、連続体を近似する有限要素法等のシミュレーション方法とを、併用することが可能である。例えば、弾性体と液体とが接触する系において、弾性体を有限要素法で解析し、液体を上記実施例による分子動力学計算により解析することが可能である。

次に、上記実施例によるシミュレーション方法を用いて三次元のダム決壊シミュレーションを行なった結果について説明する。

まず、シミュレーションの条件について説明する。使用した相互作用ポテンシャルは、レナードジョーンズ型ポテンシャルである。具体的には、下記の式で表される。

このとき、式（７）の関数ｆは、下記の式で表される。

結合定数ε、σ、ｒ_０は、下記の値とした。この値は、アルゴン（Ａｒ）原子の値に相当する。

決壊前のダムの外形を、高さＨ、横幅Ｌ、厚さＤの直方体とした。時刻ｔ＝０において、横幅方向に直交する１つの面が決壊した後の粒子系の挙動をシミュレーションした。決壊後のダムの横幅は２Ｌとした。高さＨ、横幅Ｌ、厚さＤとして、下記の値を採用した。

くりこみ変換によって粒子系の寸法が不変であるため、くりこまれた粒子系Ｓ’の決壊前のダムの外形の高さＨ’、横幅Ｌ’、厚さＤ’は、それぞれ元の粒子系Ｓの外形の高さＨ、横幅Ｌ、厚さＤと同一である。

シミュレーションの条件として、くりこみ変換後の重力加速度ｇ’、くりこみ変換後の粒子数Ｎ’、くりこみ回数ｎ、くりこみ変換後の温度Ｔ’を、下記の値とした。

重力加速度ｇ’は、くりこみ変換処理の前後において粒子系のボンド数が一致するように設定した。重力加速度ｇ’及び温度Ｔ’が上述の条件のとき粒子系Ｓ’は液体の状態である。

図６Ａ〜図６Ｄに、実施例によるシミュレーション方法を用いたシミュレーション結果を図形で示す。図６Ａ〜図６Ｄは、ダムの厚さ方向に垂直な断面を示す。図６Ａ、図６Ｂ、図６Ｃ、及び図６Ｄは、それぞれｔ＝０［ｓ］、０．００００５［ｓ］、０．０００１［ｓ］、及び０．０００２［ｓ］における粒子の分布を示す。

ダムが決壊すると、図６Ｂに示すように液体が右方に流れ出す。図６Ｃに示すように、液体が図の右側の壁に衝突すると、上方に向って壁を這い上がる。その後、図６Ｄに示すように、壁を這い上がった液体の一部が液滴となる。

図７Ａ〜図７Ｄに、特許第５２４１４６８号公報に開示された比較例によるシミュレーション方法を用いたシミュレーション結果を図形で示す。比較例においては、くりこみ変換によって粒子系の寸法が変化する。くりこまれた粒子系Ｓ’の決壊前のダムの外形の高さＨ’、横幅Ｌ’、厚さＤ’は、それぞれ下記の値になる。

シミュレーションの条件として、くりこみ変換後の重力加速度ｇ’、くりこみ変換後の粒子数Ｎ’、くりこみ回数ｎ、くりこみ変換後の温度Ｔ’を、下記の値とした。

シミュレーションされる液体のレイノルズ数及びフルード数は、実施例によるシミュレーション方法でくりこまれた粒子系Ｓ’と、比較例によるシミュレーション方法でくりこまれた粒子系Ｓ’とで同一である。このため、両者は、流体として共通した挙動を示す。

図６Ａ〜図６Ｄ、及び図７Ａ〜図７Ｄに示したように、実施例によるシミュレーション方法と、比較例によるシミュレーション方法とでは、全体としてほぼ似通った結果が得られている。ところが、細部に着目すると、両者に相違が見出される。

実施例によるシミュレーション方法では、図６Ｃに示した状態において、ケルビン−ヘルムホルツ不安定（界面の波打ち現象）が再現されている。これに対し、比較例によるシミュレーション方法では、図７Ｃに示した状態において、界面の波打ち現象が再現されていない。また、実施例によるシミュレーション方法では、図６Ｄに示した状態において、壁を這い上がった液体が落下する際に液滴が形成されている。これに対し、比較例によるシミュレーション方法では、液滴の形成が再現されていない。

上述の考察から、実施例によるシミュレーション方法を適用することにより、粒子系の挙動をより正確に再現できることがわかる。

上記実施例によるシミュレーション方法は、コンピュータがコンピュータプログラムを実行することにより実現される。コンピュータプログラムは、例えばデータ記録媒体に記録された状態で提供される。または、コンピュータプログラムは、電気通信回線を介して提供することも可能である。

図８に、実施例によるシミュレーション装置１０のブロック図を示す。シミュレーション装置１０は、入力部１１、シミュレーション処理部１２、及び作用部１３を含む。センサ２１が、対象物２０の種々の物理量、例えば、外形寸法、変形量、温度等を検出する。センサ２１の検出結果が、シミュレーション装置１０の入力部に入力される。

入力部１１は、センサ２１から入力された検出結果を、シミュレーションで使用可能なデータに変換する。シミュレーション処理部１２は、入力部１１で変換されたデータを初期条件として、図１に示したシミュレーション方法を実行する。これにより、対象物２０の将来的な挙動が推測される。

作用部１３は、シミュレーション結果に基づいて、対象物２０に対して物理的作用を施す。物理的作用がシミュレーション結果に基づいたものであるため、対象物２０に対してより適切な作用を施すことができる。

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

１０ シミュレーション装置
１１ 入力部
１２ シミュレーション処理部
１３ 作用部
２０ 対象物
２１ センサ

本発明の一観点によると、
くりこみ回数に依存するくりこみ因子αに基づいて、複数の粒子からなるシミュレーション対象の粒子系Ｓに対してくりこみ変換処理を行う工程と、
くりこまれた粒子系Ｓ’について分子動力学計算を実行することにより、くりこまれた前記粒子系Ｓ’の粒子の位置ベクトル及び運動量ベクトルを求める工程と
を有し、
前記粒子系Ｓの粒子間の相互作用ポテンシャルφが、エネルギの次元を持つ相互作用係数ε、無次元関数ｆ、粒子を特徴づけるパラメータｒ_０、σ、粒子間距離ｒを用いて、

と表すことができ、
前記粒子系Ｓの空間の次元をｄ、粒子数をＮ、各粒子の質量をｍで表し、粒子系Ｓ’の粒子数をＮ’、各粒子の質量をｍ’で表したとき、下記の変換側

を適用し、くりこまれた前記粒子系Ｓ’の相互作用ポテンシャル

に基づいて前記分子動力学計算を実行するシミュレーション方法が提供される。

分子動力学を用いたコンピュータシミュレーションが行われている。分子動力学では、シミュレーションの対象となる系を構成する粒子の運動方程式が数値的に解かれる。シミュレーション対象の系に含まれる粒子数が増加すれば、必要な計算量が増大する。現行のコンピュータの演算能力でシミュレーション可能な系の粒子数は、通常、数十万個程度である。

次に、くりこみ群の手法を応用した分子動力学（以下、くりこみ群分子動力学という。）について概念的に説明する。

くりこみ群分子動力学では、シミュレーション対象の粒子系Ｓを、粒子系Ｓの粒子数よりも少ない数の粒子で構成された粒子系Ｓ’（以下、くりこまれた粒子系Ｓ’という。）に対応づける。次に、くりこまれた粒子系Ｓ’に対して分子動力学計算を実行する。くりこまれた粒子系Ｓ’に対する計算結果を、シミュレーション対象の粒子系Ｓに対応付ける。これにより、シミュレーション対象の粒子系Ｓに対して分子動力学計算を直接的に実行する場合よりも計算量を少なくすることができる。シミュレーション対象の粒子系Ｓにおける物理量（例えば、粒子数、粒子の質量等）と、くりこまれた粒子系Ｓ’における物理量とを相互に対応付けるための変換則を、くりこみ変換則という。

積分変数ｑ_ｊについて積分を実行すれば、式（２１）を用いて下記の式が得られる。

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

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

Claims (6)

1. くりこみ回数に依存するくりこみ因子αに基づいて、複数の粒子からなるシミュレーション対象の粒子系Ｓに対してくりこみ変換処理を行う工程と、
   くりこまれた粒子系Ｓ’について分子動力学計算を実行することにより、くりこまれた前記粒子系Ｓ’の粒子の位置ベクトル及び運動量ベクトルを求める工程と
   を有し、
   前記粒子系Ｓの粒子間の相互作用ポテンシャルφが、エネルギの次元を持つ相互作用係数ε、無次元関数ｆ、粒子を特徴づけるパラメータｒ_０、σ、粒子間距離ｒを用いて、
   

   

   と表すことができ、
   前記粒子系Ｓの空間の次元をｄで表したとき、下記の変換側
   

   

   を適用し、くりこまれた前記粒子系Ｓ’の相互作用ポテンシャル
   

   

   に基づいて前記分子動力学計算を実行するシミュレーション方法。
2. 前記くりこみ変換処理を行う工程でのくりこみ回数をｎで表したとき、前記くりこみ因子αは２^ｎである請求項１に記載のシミュレーション方法。
3. 前記粒子系Ｓの温度をＴで表し、くりこまれた前記粒子系Ｓ’の温度をＴ’で表したとき、下記の変換則
   

   

   を適用して、前記分子動力学計算を実行するときの温度の初期条件を設定する請求項１または２に記載のシミュレーション方法。
4. 請求項１乃至３のいずれか１項のシミュレーション方法をコンピュータで実行するためのコンピュータプログラム。
5. 請求項４のコンピュータプログラムが、コンピュータで読み取り可能に記録された記録媒体。
6. 対象物から検出された物理量を、シミュレーションで使用可能なデータに変換する入力部と、
   前記入力部で変換された前記データを初期条件としてシミュレーションを実行するシミュレーション処理部と、
   前記シミュレーション処理部で実行されたシミュレーション結果に基づいて、前記対象物に作用を施す作用部と
   を有するシミュレーション装置であって、
   前記対象物を複数の粒子で表した粒子系Ｓの粒子間の相互作用ポテンシャルφが、エネルギの次元を持つ相互作用係数ε、無次元関数ｆ、粒子を特徴づけるパラメータｒ_０、σ、粒子間距離ｒを用いて、
   

   

   と表すことができ、
   前記シミュレーション処理部は、
   前記対象物を、複数の粒子で表した粒子系Ｓに対して、くりこみ回数に依存するくりこみ因子αに基づいてくりこみ変換処理と、
   くりこまれた粒子系Ｓ’について分子動力学計算を実行することにより、くりこまれた前記粒子系Ｓ’の粒子の位置ベクトル及び運動量ベクトルを求める処理と
   を実行し、
   前記分子動力学計算は、
   前記粒子系Ｓの空間の次元をｄで表したとき、下記の変換側
   

   

   を適用し、くりこまれた前記粒子系Ｓ’の相互作用ポテンシャル
   

   

   に基づいて実行される前記シミュレーション装置。

Images