図面 (/)

技術 シミュレーション方法、シミュレーションプログラム、及びシミュレーション装置

出願人 住友重機械工業株式会社
発明者 市嶋大路
出願日 2015年5月21日 (5年7ヶ月経過) 出願番号 2015-103371
公開日 2016年12月22日 (4年0ヶ月経過) 公開番号 2016-218767
状態 特許登録済
技術分野 特定用途計算機
主要キーワード くりこみ 相互作用係数 一次元座標 不活性原子 分配関数 単純立方格子 運動量ベクトル 三次元格子
関連する未来課題
重要な関連分野

この項目の情報は公開日時点(2016年12月22日)のものです。
また、この項目は機械的に抽出しているため、正しく解析できていない場合があります

図面 (9)

課題

粒子速度分布が維持されるくりこみ変換を適用してシミュレーションを行う方法を提供する。

解決手段

くりこみ回数に依存するくりこみ因子αに基づいて、シミュレーション対象粒子系Sに対してくりこみ変換処理を行う。くりこまれた粒子系S’について分子動力学計算を実行することにより、粒子系S’の粒子の位置ベクトル及び運動量ベクトルを求める。粒子系Sの粒子間の相互作用ポテンシャルφが、エネルギ次元を持つ相互作用係数ε無次元関数f、粒子を特徴づけるパラメータr0、σ、粒子間距離rを用いて、φ(r)=εf((r−r0)/σ)と表される。粒子系Sの空間の次元をdで表したとき、変換側N’=N/αd、m’=mαd、ε’=εαd、ro’=αr0、σ’=ασを適用し、粒子系S’の相互作用ポテンシャルφ’(r)=ε’f((r−r0’)/σ’)に基づいて分子動力学計算を実行する。

概要

背景

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

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

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

第1のくりこみ因子α、第2のくりこみ因子γ、及び第3のくりこみ因子δが決定される。第1のくりこみ因子αは1より大きい。第2のくりこみ因子γは、0以上であり、かつ空間の次元数d以下である。第3のくりこみ因子δは、0以上である。くりこみの回数をnで表したとき、第1のくりこみ因子αは、
α=2n
と表される。

くりこみ手法を用いてくりこみ変換された粒子系S’の粒子数をN’、各粒子の質量をm’、相互作用係数をε’で表す。以下の変換式を用いて、くりこみ変換された粒子系S’の粒子数N’、質量m’、及び相互作用係数ε’を求める。
N’=N/αd
m’=mαδ−γ
ε’=εαγ

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

概要

粒子の速度分布が維持されるくりこみ変換を適用してシミュレーションを行う方法を提供する。くりこみ回数に依存するくりこみ因子αに基づいて、シミュレーション対象の粒子系Sに対してくりこみ変換処理を行う。くりこまれた粒子系S’について分子動力学計算を実行することにより、粒子系S’の粒子の位置ベクトル及び運動量ベクトルを求める。粒子系Sの粒子間の相互作用ポテンシャルφが、エネルギの次元を持つ相互作用係数ε、無次元関数f、粒子を特徴づけるパラメータr0、σ、粒子間距離rを用いて、φ(r)=εf((r−r0)/σ)と表される。粒子系Sの空間の次元をdで表したとき、変換側N’=N/αd、m’=mαd、ε’=εαd、ro’=αr0、σ’=ασを適用し、粒子系S’の相互作用ポテンシャルφ’(r)=ε’f((r−r0’)/σ’)に基づいて分子動力学計算を実行する。

目的

本発明の目的は、粒子の速度分布が維持されるくりこみ変換を適用してシミュレーションを行う方法を提供する

効果

実績

技術文献被引用数
0件
牽制数
0件

この技術が所属する分野

(分野番号表示ON)※整理標準化データをもとに当社作成

ライセンス契約や譲渡などの可能性がある特許掲載中! 開放特許随時追加・更新中 詳しくはこちら

請求項1

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

請求項2

前記くりこみ変換処理を行う工程でのくりこみ回数をnで表したとき、前記くりこみ因子αは2nである請求項1に記載のシミュレーション方法。

請求項3

前記粒子系Sの温度をTで表し、くりこまれた前記粒子系S’の温度をT’で表したとき、下記の変換則を適用して、前記分子動力学計算を実行するときの温度の初期条件を設定する請求項1または2に記載のシミュレーション方法。

請求項4

請求項1乃至3のいずれか1項のシミュレーション方法をコンピュータで実行するためのコンピュータプログラム

請求項5

請求項4のコンピュータプログラムが、コンピュータで読み取り可能に記録された記録媒体

請求項6

対象物から検出された物理量を、シミュレーション使用可能なデータに変換する入力部と、前記入力部で変換された前記データを初期条件としてシミュレーションを実行するシミュレーション処理部と、前記シミュレーション処理部で実行されたシミュレーション結果に基づいて、前記対象物に作用を施す作用部とを有するシミュレーション装置であって、前記対象物を複数の粒子で表した粒子系Sの粒子間の相互作用ポテンシャルφが、エネルギの次元を持つ相互作用係数ε、無次元関数f、粒子を特徴づけるパラメータr0、σ、粒子間距離rを用いて、と表すことができ、前記シミュレーション処理部は、前記対象物を、複数の粒子で表した粒子系Sに対して、くりこみ回数に依存するくりこみ因子αに基づいてくりこみ変換処理と、くりこまれた粒子系S’について分子動力学計算を実行することにより、くりこまれた前記粒子系S’の粒子の位置ベクトル及び運動量ベクトルを求める処理とを実行し、前記分子動力学計算は、前記粒子系Sの空間の次元をdで表したとき、下記の変換側を適用し、くりこまれた前記粒子系S’の相互作用ポテンシャルに基づいて実行される前記シミュレーション装置。

技術分野

背景技術

0002

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

0003

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

0004

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

0005

第1のくりこみ因子α、第2のくりこみ因子γ、及び第3のくりこみ因子δが決定される。第1のくりこみ因子αは1より大きい。第2のくりこみ因子γは、0以上であり、かつ空間の次元数d以下である。第3のくりこみ因子δは、0以上である。くりこみの回数をnで表したとき、第1のくりこみ因子αは、
α=2n
と表される。

0006

くりこみ手法を用いてくりこみ変換された粒子系S’の粒子数をN’、各粒子の質量をm’、相互作用係数をε’で表す。以下の変換式を用いて、くりこみ変換された粒子系S’の粒子数N’、質量m’、及び相互作用係数ε’を求める。
N’=N/αd
m’=mαδ−γ
ε’=εαγ

0007

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

先行技術

0008

特許第5241468号公報

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

0009

くりこみ変換された粒子系S’のマクスウェル速度分布則fmax(v’)は、以下の式で表される。
fmax(v’)=exp[(−m’/2kBT)v’2]
ここで、v’は、くりこみ変換された粒子系S’の各粒子の速度であり、kBはボルツマン定数であり、Tは粒子系S’の温度である。

0010

上述のマクスウェル速度分布則fmax(v’)は、元の粒子系Sの各粒子の質量mを用いて、以下の式に書き直される。
fmax(v’)=exp[(−m/2kBT)(v’α(δ−γ)/2)2]
上述の式は、くりこみ変換された粒子系S’の粒子の速度分布の標準偏差が、元の粒子系Sの粒子の速度分布の標準偏差の1/α(δ−γ)/2倍になっていることを意味する。すなわち、くりこみ変換された粒子系S’では、元の粒子系Sの粒子の速度分布が再現されていない。くりこみ回数nが大きくなると、くりこみ変換された粒子系S’の粒子の速度分布が、元の粒子系Sの粒子の速度分布から大きく乖離してしまう。

0011

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

0012

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

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

0013

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



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



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



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

0014

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

発明の効果

0015

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

図面の簡単な説明

0016

図1は、実施例によるシミュレーション方法のフローチャートである。
図2は、一次元状に配置された粒子を示す概念図である。
図3は、相互作用ポテンシャルφがレナードジョーンズポテンシャルである場合の式(21)の算出結果、及び式(22)の数値積分結果を示すグラフである。
図4は、相互作用ポテンシャルφがモース型ポテンシャルである場合の式(21)の算出結果、及び式(22)の数値積分結果を示すグラフである。
図5A〜図5Cは、二次元格子状に配置された粒子、及び粒子間相互作用を示す概念図である。
図6A〜図6Dは、実施例によるシミュレーション方法を用いたシミュレーション結果を示す図である。
図7A〜図7Dは、比較例によるシミュレーション方法を用いたシミュレーション結果を示す図である。
図8は、実施例によるシミュレーション装置のブロック図である。

実施例

0017

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



ここで、
mは粒子の質量を表し、
φは粒子間の相互作用ポテンシャルを表し、
ベクトルpjは粒子の運動量ベクトルを表し、
ベクトルqjは粒子の位置ベクトル(位置座標)を表す。

0018

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

0019

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

0020

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

0021

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

0022

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



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

0023

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

0024

式(8)は、下記の式のように変形することができる。



式(9)からわかるように、レナードジョーンズ型ポテンシャルは(r−r0)/σの関数であり、式(7)の形に表すことができる。

0025

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

0026

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



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

0027

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

0028

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

0029

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

0030

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

0031

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

0032

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

0033

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

0034

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

0035

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

0036

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



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



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

0037

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

0038

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

0039

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

0040

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



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

0041

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

0042

積分変数qjについて積分を実行すれば、式(21)を用いて下記の式が得られる。



ここで、raは粒子の直径を表し、z(qi−qk)及びP(qi−qk)は、下記の式で表される。






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

0043

次に、z(qi−qk)を具体的に計算する。相互作用ポテンシャルφがレナードジョーンズ型ポテンシャルの場合、φ(2n)は、下記の式で表される。



相互作用ポテンシャルφがモース型ポテンシャルの場合、φ(2n)は、下記の式で表される。

0044

式(25)または式(26)を式(24)に代入して数値積分を行う。式(24)に式(25)または式(26)を代入する際には、式(20)を用いる。数値積分において、式(24)の積分範囲に現れているaは、近似的にr0に等しいと仮定する。

0045

図3に、相互作用ポテンシャルφがレナードジョーンズ型ポテンシャルである場合の式(23)の算出結果、及び式(24)の数値積分結果を示す。図4に、相互作用ポテンシャルφがモース型ポテンシャルである場合の式(23)の算出結果、及び式(24)の数値積分結果を示す。粒子i、粒子j、粒子k、粒子l、及び粒子mがこの順番に一次元状に配置されてい場合の粒子kの位置座標がqkで表される。図3及び図4横軸は、qk/2を表し、縦軸は、P(qi−qk)P(qk−qm)、及びz(qi−qk)z(qk−qm)を、対数目盛で表す。図3に至る数値計算において、ε/kBT=2.0、r0/σ=1.12とした。図4に至る数値計算において、ε/kBT=2.0、r0/σ=2.24とした。

0046

相互作用ポテンシャルφがレナードジョーンズ型ポテンシャルである場合、及びモース型ポテンシャルである場合のいずれにおいても、P(qi−qk)P(qk−qm)の変化に比べて、z(qi−qk)z(qk−qm)の変化が緩やかであることがわかる。このため、P(qi−qk)P(qk−qm)に対してz(qi−qk)z(qk−qm)は、ほぼ定数であると近似できる。

0047

粒子kが位置座標qkに存在する確率p(qk)は、以下のように近似することができる。

0048

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

0049

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

0050

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

0051

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

0052

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

0053

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

0054

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



ここで、<i,j>は、最近接格子間で和を取ることを意味する。

0055

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

0056

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



式(31)の導出に際し、下記の式が利用される。ここで、運動量ベクトルpj2は、ベクトルの内積を意味する。

0057

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

0058

式(18)に、式(30)及び式(31)を代入することにより、結果に影響しない係数を除いて下記の式が得られる。

0059

式(33)から、粗視化されたハミルトニアン(くりこまれた粒子系S’のハミルトニアン)H’は、以下の式で表される。

0060

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

0061

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

0062

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

0063

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



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

0064

式(38)から、式(11)に示したくりこみ変換測、及び式(14)に示したくりこまれた粒子系S’のハミルトニアンH’が導き出される。

0065

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

0066

式(40)に、式(11)のm’に関する式、及び式(15)を代入することにより、下記の式が得られる。

0067

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

0068

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

0069

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

0070

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

0071

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



このとき、式(7)の関数fは、下記の式で表される。

0072

結合定数ε、σ、r0は、下記の値とした。この値は、アルゴン(Ar)原子の値に相当する。

0073

決壊前のダムの外形を、高さH、横幅L、厚さDの直方体とした。時刻t=0において、横幅方向に直交する1つの面が決壊した後の粒子系の挙動をシミュレーションした。決壊後のダムの横幅は2Lとした。高さH、横幅L、厚さDとして、下記の値を採用した。



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

0074

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



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

0075

図6A〜図6Dに、実施例によるシミュレーション方法を用いたシミュレーション結果を図形で示す。図6A〜図6Dは、ダムの厚さ方向に垂直な断面を示す。図6A、図6B、図6C、及び図6Dは、それぞれt=0[s]、0.00005[s]、0.0001[s]、及び0.0002[s]における粒子の分布を示す。

0076

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

0077

図7A〜図7Dに、特許第5241468号公報に開示された比較例によるシミュレーション方法を用いたシミュレーション結果を図形で示す。比較例においては、くりこみ変換によって粒子系の寸法が変化する。くりこまれた粒子系S’の決壊前のダムの外形の高さH’、横幅L’、厚さD’は、それぞれ下記の値になる。

0078

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

0079

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

0080

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

0081

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

0082

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

0083

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

0084

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

0085

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

0086

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

0087

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

0088

10シミュレーション装置
11 入力部
12シミュレーション処理部
13 作用部
20対象物
21 センサ

ページトップへ

この技術を出願した法人

この技術を発明した人物

ページトップへ

関連する挑戦したい社会課題

関連する公募課題

ページトップへ

技術視点だけで見ていませんか?

この技術の活用可能性がある分野

分野別動向を把握したい方- 事業化視点で見る -

(分野番号表示ON)※整理標準化データをもとに当社作成

ページトップへ

おススメ サービス

おススメ astavisionコンテンツ

新着 最近 公開された関連が強い技術

この 技術と関連性が強い人物

関連性が強い人物一覧

この 技術と関連する社会課題

関連する挑戦したい社会課題一覧

この 技術と関連する公募課題

関連する公募課題一覧

astavision 新着記事

サイト情報について

本サービスは、国が公開している情報(公開特許公報、特許整理標準化データ等)を元に構成されています。出典元のデータには一部間違いやノイズがあり、情報の正確さについては保証致しかねます。また一時的に、各データの収録範囲や更新周期によって、一部の情報が正しく表示されないことがございます。当サイトの情報を元にした諸問題、不利益等について当方は何ら責任を負いかねることを予めご承知おきのほど宜しくお願い申し上げます。

主たる情報の出典

特許情報…特許整理標準化データ(XML編)、公開特許公報、特許公報、審決公報、Patent Map Guidance System データ