(57)【要約】【課題】例えば気体状態のように粒子密度が低く、かつ流動のある粒子系の解析においても、温度制御を行うことが可能なシミュレーション方法を提供する。【解決手段】解析領域に複数の粒子を配置して、複数の前記粒子からなる粒子系の挙動を、分子動力学法を用いてシミュレーションする際に、粒子ごとに、着目する粒子を含む局所領域の内側に存在する粒子の重心速度を、局所領域の内側に存在する粒子の数に基づいて変換して局所領域の重心速度の変換値を算出する。粒子の速度と、当該粒子を含む局所領域の重心速度の変換値との差を、粒子系の温度を算出する基礎となる粒子の速度として用い、粒子系の温度を目標温度に維持する温度制御を行う。【選択図】図３
2019070896.tif 000002
【特許請求の範囲】
【請求項１】
解析領域に複数の粒子を配置して、複数の前記粒子からなる粒子系の挙動を、分子動力学法を用いてシミュレーションする方法であって、
前記粒子ごとに、着目する粒子を含む局所領域の内側に存在する前記粒子の重心速度を、前記局所領域の内側に存在する前記粒子の数に基づいて変換して前記局所領域の重心速度の変換値を算出し、
前記粒子の速度と、当該粒子を含む前記局所領域の重心速度の変換値との差を、前記粒子系の温度を算出する基礎となる前記粒子の速度として用い、前記粒子系の温度を目標温度に維持する温度制御を行うシミュレーション方法。
【請求項２】
前記局所領域の重心速度の変換値は、前記局所領域の内側に存在する前記粒子の数が増加するに従って単調に増加し、前記局所領域の内側に存在する前記粒子の重心速度に近づく請求項１に記載のシミュレーション方法。
【請求項３】
前記局所領域の内側に存在する前記粒子の数が１のとき前記局所領域の重心速度の変換値が０であり、前記局所領域の内側に存在する前記粒子の数が増加するに従って前記局所領域の内側に存在する前記粒子の重心速度に漸近する請求項１または２に記載のシミュレーション方法。
【請求項４】
複数の前記粒子のうち着目する粒子ｉを前記局所領域の内側に存在する前記粒子の重心速度をＶｇ，ｉで表し、前記局所領域の内側に存在する前記粒子の数をｎｉで表したとき、重心速度の変換値Ｗｇ，ｉを、
Ｗｇ，ｉ＝（１−（１／ｎｉ２））Ｖｇ，ｉ
と定義する請求項１乃至３のいずれか１項に記載のシミュレーション方法。
【請求項５】
解析領域に複数の粒子を配置して、複数の前記粒子からなる粒子系の挙動を、分子動力学法を用いてシミュレーションする処理手段を有するシミュレーション装置であって、
前記処理手段は、
入力手段から、シミュレーションの初期条件、境界条件、及び前記粒子系の目標温度を取得する機能、
入力された前記初期条件、前記境界条件に基づいて、前記粒子系の挙動をシミュレーションする機能、
複数の前記粒子の各々について、着目する粒子を含む局所領域の内側に存在する前記粒子の重心速度を、前記局所領域の内側に存在する前記粒子の数に基づいて変換して前記局所領域の重心速度の変換値を算出し、前記粒子の速度と、当該粒子を含む前記局所領域の重心速度の変換値との差を、前記粒子系の温度を算出する基礎となる前記粒子の速度として、前記粒子系の温度が前記目標温度に維持されるように温度制御を行う機能、及び
シミュレーション結果を出力手段に出力する機能
を有するシミュレーション装置。
【請求項６】
前記局所領域の重心速度の変換値は、前記局所領域の内側に存在する前記粒子の数が増加するに従って単調に増加し、前記局所領域の内側に存在する前記粒子の重心速度に近づく請求項５に記載のシミュレーション装置。
【請求項７】
前記局所領域の内側に存在する前記粒子の数が１のとき前記局所領域の重心速度の変換値が０であり、前記局所領域の内側に存在する前記粒子の数が増加するに従って前記局所領域の内側に存在する前記粒子の重心速度に漸近する請求項５または６に記載のシミュレーション装置。
【請求項８】
複数の前記粒子のうち着目する粒子ｉを前記局所領域の内側に存在する前記粒子の重心速度をＶｇ，ｉで表し、前記局所領域の内側に存在する前記粒子の数をｎｉで表したとき、重心速度の変換値Ｗｇ，ｉを、
Ｗｇ，ｉ＝（１−（１／ｎｉ２））Ｖｇ，ｉ
と定義する請求項５乃至７のいずれか１項に記載のシミュレーション装置。
【請求項９】
解析領域に複数の粒子を配置して、複数の前記粒子からなる粒子系の挙動を、分子動力学法を用いてシミュレーションする機能をコンピュータに実現させるプログラムであって、
与えられた初期条件、境界条件に基づいて、前記粒子系の挙動をシミュレーションする機能、及び
複数の前記粒子の各々について、着目する粒子を含む局所領域の内側に存在する前記粒子の重心速度を、前記局所領域の内側に存在する前記粒子の数に基づいて変換して前記局所領域の重心速度の変換値を算出し、前記粒子の速度と、当該粒子を含む前記局所領域の重心速度の変換値との差を、前記粒子系の温度を算出する基礎となる前記粒子の速度として、前記粒子系の温度が目標温度に維持されるように温度制御を行う機能
を有するプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、シミュレーション方法、シミュレーション装置、及びプログラムに関する。
【背景技術】
【０００２】
分子動力学法に温度制御を取り入れてシミュレーションを行う方法が公知である（例えば、特許文献１）。例えば、分子動力学法によって算出された粒子の速度を、所望の設定温度、及び系の温度に基づいて更新する。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特開平１０−３３４０７６号公報
【発明の概要】
【発明が解決しようとする課題】
【０００４】
分子運動論において、複数の分子からなる系の温度と、複数の分子の運動エネルギの平均値とが、ボルツマン定数によって関連付けられる。温度は慣性系に依存しないことから、流動のある系の解析を行う場合には、系の温度を算出するにあたり、各分子の速度から系の重心速度を除いた値に基づいて、分子の運動エネルギを算出する必要がある。
【０００５】
一般的に、ある粒子の位置における重心速度を求める場合、当該粒子の周辺（局所領域）の粒子群の平均速度から重心速度を求める。ところが、複数の粒子が気体状態になっている場合、粒子密度が低いため、１つの粒子の周辺に他の粒子が存在しなくなる。すなわち、局所領域内に、着目する１つの粒子しか存在しなくなる。この場合、着目する粒子の位置における重心速度が、その粒子の速度と等しくなる。粒子の速度と、その位置における重心速度との差がゼロになるため、算出される温度もゼロになってしまう。このため、適切な温度制御を行うことができない。
【０００６】
本発明の目的は、例えば気体状態のように粒子密度が低く、かつ流動のある粒子系の解析においても、温度制御を行うことが可能なシミュレーション方法、シミュレーション装置、及びプログラムを提供することである。
【課題を解決するための手段】
【０００７】
本発明の一観点によると、
解析領域に複数の粒子を配置して、複数の前記粒子からなる粒子系の挙動を、分子動力学法を用いてシミュレーションする方法であって、
前記粒子ごとに、着目する粒子を含む局所領域の内側に存在する前記粒子の重心速度を、前記局所領域の内側に存在する前記粒子の数に基づいて変換して前記局所領域の重心速度の変換値を算出し、
前記粒子の速度と、当該粒子を含む前記局所領域の重心速度の変換値との差を、前記粒子系の温度を算出する基礎となる前記粒子の速度として用い、前記粒子系の温度を目標温度に維持する温度制御を行うシミュレーション方法が提供される。
【０００８】
本発明の他の観点によると、
解析領域に複数の粒子を配置して、複数の前記粒子からなる粒子系の挙動を、分子動力学法を用いてシミュレーションする処理手段を有するシミュレーション装置であって、
前記処理手段は、
入力手段から、シミュレーションの初期条件、境界条件、及び前記粒子系の目標温度を取得する機能、
入力された前記初期条件、前記境界条件に基づいて、前記粒子系の挙動をシミュレーションする機能、
複数の前記粒子の各々について、着目する粒子を含む局所領域の内側に存在する前記粒子の重心速度を、前記局所領域の内側に存在する前記粒子の数に基づいて変換して前記局所領域の重心速度の変換値を算出し、前記粒子の速度と、当該粒子を含む前記局所領域の重心速度の変換値との差を、前記粒子系の温度を算出する基礎となる前記粒子の速度として、前記粒子系の温度が前記目標温度に維持されるように温度制御を行う機能、及び
シミュレーション結果を出力手段に出力する機能
を有するシミュレーション装置が提供される。
【０００９】
本発明のさらに他の観点によると、
解析領域に複数の粒子を配置して、複数の前記粒子からなる粒子系の挙動を、分子動力学法を用いてシミュレーションする機能をコンピュータに実現させるプログラムであって、
与えられた初期条件、境界条件に基づいて、前記粒子系の挙動をシミュレーションする機能、及び
複数の前記粒子の各々について、着目する粒子を含む局所領域の内側に存在する前記粒子の重心速度を、前記局所領域の内側に存在する前記粒子の数に基づいて変換して前記局所領域の重心速度の変換値を算出し、前記粒子の速度と、当該粒子を含む前記局所領域の重心速度の変換値との差を、前記粒子系の温度を算出する基礎となる前記粒子の速度として、前記粒子系の温度が目標温度に維持されるように温度制御を行う機能
を有するプログラムが提供される。
【発明の効果】
【００１０】
例えば気体状態のように粒子密度が低く、かつ流動のある粒子系の解析においても、適切に温度制御を行うことが可能である。
【図面の簡単な説明】
【００１１】
【図１】図１Ａは、粒子ｉ、及びその周辺（粒子ｉの位置を中心とする局所領域Ｓｉ）に位置する複数の粒子ｊの速度を示す模式図であり、図１Ｂは、局所領域Ｓｉの内側に存在する粒子ｊの、重心速度Ｖｇ，ｉに対する相対速度ｖｊ−Ｖｇ，ｉを示す模式図である。
【図２】図２Ａは、気体状態のように粒子密度が低い状態の場合の粒子ｉ、及び局所領域Ｓｉの模式図であり、図２Ｂは、局所領域Ｓｉの内側に存在する粒子ｉの、重心速度Ｖｇ，ｉに対する相対速度ｖｉ−Ｖｇ，ｉを示す模式図である。
【図３】図３は、実施例によるシミュレーション方法のフローチャートである。
【図４】図４は、シミュレーション対象の粒子系の模式図である。
【図５】図５は、シミュレーションによって算出された粒子系の温度の時間変化を、比較例と比較して示すグラフである。
【図６】図６は、実施例によるシミュレーション装置のブロック図である。
【発明を実施するための形態】
【００１２】
実施例について説明する前に、図１Ａ〜図２Ｂを参照して、分子動力学法のシミュレーション対象である複数の粒子からなる粒子系の一般的な温度の算出方法、及びその課題について説明する。
【００１３】
分子運動論において、分子の運動エネルギと温度との関係は、以下の式で表される。
【００１４】
【数１】
2019070896.tif 000003
【００１５】
ここで、Ｎは粒子数、ｍｉは粒子ｉの質量、ｖｉは粒子ｉの速度、Ｔは粒子系の温度、ｋＢはボルツマン定数である。分子動力学法では、目標とする温度の粒子系を実現するために、式（１）を用いて粒子系の温度Ｔを求めることが多い。
【００１６】
系の温度は慣性系に依存しないため、流動等の、重心が移動する粒子系の解析を行う場合には、粒子系の温度Ｔを求める際に、各粒子の速度から粒子系の重心速度を除いた値を用いる必要がある。
【００１７】
図１Ａは、粒子ｉ、及びその周辺（粒子ｉの位置を中心とする局所領域Ｓｉ）に位置する複数の粒子ｊの速度を示す模式図である。局所領域Ｓｉの内側に存在する粒子ｊ（粒子ｉを含む）の個数をｎｉで表し、粒子ｊの速度をｖｊで表すと、局所領域Ｓｉの内側に存在する粒子ｊの重心速度（以下、局所領域Ｓｉの重心速度という。）Ｖｇ，ｉは、以下の式で計算することができる。
【００１８】
【数２】
2019070896.tif 000004
【００１９】
図１Ｂは、局所領域Ｓｉの内側に存在する粒子ｊの、重心速度Ｖｇ，ｉに対する相対速度ｖｊ−Ｖｇ，ｉを示す模式図である。相対速度ｖｊ−Ｖｇ，ｉの平均はゼロである。粒子系の温度Ｔは、以下の式に基づいて求めることができる。
【００２０】
【数３】
2019070896.tif 000005
【００２１】
式（３）は、粒子系の温度Ｔを求める際に、粒子ｉの速度として、粒子ｉの速度ｖｉから粒子ｉの位置を中心とする局所領域Ｓｉの重心速度Ｖｇ，ｉを差し引いた値を使用することを意味している。
【００２２】
図２Ａは、気体状態のように粒子密度が低い状態の場合の粒子ｉ、及び局所領域Ｓｉの模式図である。粒子密度が低いため、局所領域Ｓｉの内側に粒子ｉ以外の粒子が存在しない。この時の局所領域Ｓｉの重心速度Ｖｇ，ｉは、粒子ｉの速度ｖｉに等しい。
【００２３】
図２Ｂは、局所領域Ｓｉの内側に存在する粒子ｉの、重心速度Ｖｇ，ｉに対する相対速度ｖｉ−Ｖｇ，ｉを示す模式図である。相対速度ｖｉ−Ｖｇ，ｉはゼロになる。このときの粒子系の温度Ｔを式（３）から求めると、温度Ｔがゼロになってしまう。局所領域Ｓｉの内側の粒子数ｎｉが１ではない場合でも、粒子数ｎｉが少ない場合には、温度Ｔが実際の温度より低く見積もられてしまう。このため、分子動力学法によるシミュレーションにおいて、適切な温度制御を行うことができない。
【００２４】
次に、図３〜図５を参照して、適切な温度制御を行うことが可能な実施例によるシミュレーション方法について説明する。
【００２５】
図３は、実施例によるシミュレーション方法のフローチャートである。まず、シミュレーションの初期条件、境界条件、及び目標温度を決定する（ステップＳ１）。初期条件には、粒子系の各粒子の位置及び速度が含まれる。境界条件として、例えば周期境界条件を適用する。シミュレーションにおいて、粒子系の温度Ｔが目標温度に近づくように温度制御される。
【００２６】
次に、粒子が従う運動方程式に基づいて動解析を行う（ステップＳ２）。粒子間の相互作用ポテンシャルとして、例えばレナードジョーンズポテンシャルを用いることができる。この動解析により、１タイムステップ経過後の粒子の位置及び速度が求まる。
【００２７】
次に、粒子系を構成する粒子ごとに、着目する粒子ｉを含む局所領域Ｓｉの重心速度Ｖｇ，ｉを、局所領域Ｓｉの内側に存在する粒子の数ｎｉに基づいて変換して局所領域Ｓｉの重心速度の変換値Ｗｇ，ｉを算出する（ステップＳ３）。例えば、重心速度の変換値Ｗｇ，ｉを以下の変換則に基づいて算出する。
【００２８】
【数４】
2019070896.tif 000006
【００２９】
式（４）の変換則を用いると、局所領域Ｓｉの内側に存在する粒子数ｎｉが少なくなるに従って、重心速度の変換値Ｗｇ，ｉがゼロに近づく。粒子数ｎｉが１のとき、重心速度の変換値Ｗｇ，ｉはゼロになる。粒子数ｎｉが大きくなると、重心速度の変換値Ｗｇ，ｉは重心速度Ｖｇ，ｉに近づく。粒子数ｎｉが無限大のとき、重心速度の変換値Ｗｇ，ｉは重心速度Ｖｇ，ｉに等しくなる。
【００３０】
次に、局所領域Ｓｉの重心速度の変換値Ｗｇ，ｉを用いて温度制御を行う（ステップＳ４）。粒子系の温度Ｔは、以下の式で算出される。
【００３１】
【数５】
2019070896.tif 000007
【００３２】
温度制御には、例えばランジュバン法を用いることができる。ランジュバン法においては、粒子にランダム力と散逸力とを与えて温度制御を行う。ランジュバン法については、Ｊ．Ｍ．テッセン著、「計算物理学」丸善出版（２００３年）の２０４〜２０８頁に説明されている。
【００３３】
実行すべきタイムステップ数だけ、ステップＳ２からステップＳ４まで繰り返した後、シミュレーション結果を出力する（ステップＳ５）。実行すべきタイムステップ数は、予め設定されている。
【００３４】
実際に、水滴と水蒸気とが混在する粒子系についてシミュレーションを行った。
【００３５】
図４は、シミュレーション対象の粒子系の模式図である。立方体の解析領域１０内に複数の粒子１１が配置されている。解析領域１０のほぼ中心に、粒子密度の高い領域が存在する。粒子密度の高い領域は、水滴１２に対応する。水滴１２の周囲の粒子密度の低い領域には、水蒸気状態の粒子１１が存在する。解析領域１０の前後、上下、左右の境界には、周期境界条件を適用した。粒子数を８０００個とし、粒子系の目標温度を３００Ｋとした。
【００３６】
粒子間の相互作用ポテンシャルＵ（ｒ）として、下記のレナードジョーンズポテンシャルを用いた。
【００３７】
【数６】
2019070896.tif 000008
【００３８】
ここで、ε及びσはフィッティングパラメータである。本シミュレーションでは、ε及びσとして水分子を想定した値を用いた。具体的には、εを４０４．５２Ｋとし、σを０．２６４ｎｍとした。式（５）の質量ｍｉを１８．０２Ｄａとした。抵抗係数は０．１τとした。ここで、τは以下の式で定義される。
【００３９】
【数７】
2019070896.tif 000009
【００４０】
局所領域Ｓｉとして、粒子が実質的に相互作用を及ぼす距離を半径とする球体を採用した。
【００４１】
図５は、シミュレーションによって算出された粒子系の温度Ｔの時間変化を、比較例と比較して示すグラフである。横軸は経過時間を単位「ｐｓ」で表し、縦軸は温度Ｔを単位「Ｋ」で表す。図５のグラフにおいて、実線は実施例による方法でシミュレーションした結果を示し、破線は比較例による方法でシミュレーションした結果を示す。比較例においては、重心速度Ｖｇ，ｉを変換することなく式（３）を用いて粒子系の温度Ｔを算出した。
【００４２】
比較例においては、目標温度を３００Ｋとして粒子系の温度Ｔを目標温度に維持する温度制御を行っているにもかかわらず、温度Ｔが時間の経過とともに上昇し、目標温度から離れてしまうことが確認される。これは、図２Ａ及び図２Ｂを参照して説明したように、粒子系の温度Ｔが実際の温度より低く見積もられてしまうためである。
【００４３】
本実施例によるシミュレーションを用いると、粒子系の温度Ｔがほぼ目標温度３００Ｋに維持されていることがわかる。このように、実施例によるシミュレーションを用いることにより、適切な温度制御を行うことができる。
【００４４】
次に、上記実施例の変形例について説明する。上記実施例では、粒子系の温度Ｔを算出する基礎となる重心速度の変換値Ｗｇ，ｉとして式（４）を用いたが、その他の変換則により重心速度の変換値Ｗｇ，ｉを算出してもよい。例えば、局所領域Ｓｉの重心速度の変換値Ｗｇ，ｉが、局所領域Ｓｉの内側に存在する粒子の数ｎｉが増加するに従って単調に増加し、局所領域Ｓｉの重心速度Ｖｇ，ｉに近づくような変換則を用いるとよい。さらに、局所領域Ｓｉの内側に存在する粒子の数ｎｉが１のとき重心速度の変換値Ｗｇ，ｉが０になり、粒子の数ｎｉが増加するに従って局所領域Ｓｉの重心速度Ｖｇ，ｉに漸近するような変換則を用いるとよい。式（３）の変換則は、このような変換則の一例である。
【００４５】
上記実施例では、温度制御にランジュバン法を用いたが、他の方法を用いてもよい。例えば、粒子系の温度Ｔが目標温度に維持されるように粒子の速度をスケーリングしてもよい。その他に、物理的に厳密な熱浴を与える能勢フーバー法を用いてもよい。
【００４６】
上記実施例では、局所領域Ｓｉとして、粒子ｉが実質的に相互作用を及ぼす距離を半径とする球体を採用したが、その他の範囲を局所領域Ｓｉとして採用してもよい。例えば、ステップＳ２（図３）の動解析を行うときに、着目する粒子ｉと相互作用を及ぼす範囲として採用される範囲を局所領域Ｓｉとするとよい。なお、局所領域Ｓｉは、粒子ｉが実質的に他の粒子と相互作用を及ぼす範囲を含むようにするとよい。
【００４７】
分子動力学法をマクロスケールの系を扱えるように発展させた繰り込み群分子動力学法（ＲＭＤ法）が提案されている。上記実施例によるシミュレーション方法は、繰り込み群分子動力学法を用いたシミュレーションにも適用可能である。繰り込み群分子動力学法は、分子動力学法の一形態を考えることができ、本明細書において「分子動力学法」は、「繰り込み群分子動力学法」を含む上位概念として用いられている。
【００４８】
次に、図６を参照して、実施例によるシミュレーション装置について説明する。実施例によるシミュレーション装置は、図３に示したシミュレーション方法をコンピュータに実行させるプログラムを含み、このプログラムを実行する。
【００４９】
図６は、実施例によるシミュレーション装置のブロック図である。このシミュレーション装置として、シミュレーションプログラムを実行するコンピュータを用いることができる。このコンピュータは、中央処理装置（ＣＰＵ）２０、メインメモリ２１、補助記憶装置２２、入力装置２３、及び出力装置２４を含む。補助記憶装置２２は、シミュレーションプログラムが格納されている記録媒体を含む。この記録媒体は、補助記憶装置２２に内蔵されたものでもよく、補助記憶装置２２に対して着脱可能なリムーバブル媒体でもよい。このシミュレーションプログラムがメインメモリ２１にロードされ、ＣＰＵ２０によって実行される。
【００５０】
入力装置２３からシミュレーションに必要な情報が入力される。例えば、ステップＳ１（図３）で設定される種々の情報が入力される。ＣＰＵ２０は、入力装置２３から入力された初期条件、境界条件、及び目標温度を取得する。さらに、ＣＰＵ２０は、シミュレーション結果を出力装置２４に出力する。例えば、ステップＳ５（図３）において、粒子の挙動を可視化した画像、図５に示した粒子系の温度の時間変化を示すグラフ等を出力装置２４に表示する。
【００５１】
上記実施例及び変形例は例示であり、実施例及び変形例で示した構成の部分的な置換または組み合わせが可能であることは言うまでもない。実施例及び変形例の同様の構成による同様の作用効果については実施例及び変形例ごとには逐次言及しない。さらに、本発明は上述の実施例及び変形例に制限されるものではない。例えば、種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。
【符号の説明】
【００５２】
１０  解析領域
１１  粒子
１２  水滴
２０  中央処理装置（ＣＰＵ）
２１  メインメモリ
２２  補助記憶装置
２３  入力装置
２４  出力装置
【図１】
2019070896.tif 000010
【図２】
2019070896.tif 000011
【図３】
2019070896.tif 000012
【図４】
2019070896.tif 000013
【図５】
2019070896.tif 000014
【図６】
2019070896.tif 000015