(57)【特許請求の範囲】
【請求項１】
流れ場を形成する複数の粒子を含む粒子系を分子動力学法によりシミュレーションする方法であって、
着目している粒子の速度と、着目している粒子の位置における流れ場の流速との差に基づく運動エネルギ、及び着目している粒子の周囲に存在する他の複数の粒子の速度と、前記流速との差に基づく運動エネルギの平均を算出し、
算出された平均値に基づいて、粒子の各々の位置における温度を算出し、
粒子ごとに、各粒子の位置における算出された温度と、前記粒子系の目標温度とに基づいて、温度制御を行うことにより、前記粒子系の温度を目標温度に近づけるシミュレーション方法。
【請求項２】
粒子の位置における流れ場の流速は、着目している粒子の速度、及び着目している粒子の周囲に存在する他の複数の粒子の速度を平均することにより求める請求項１に記載のシミュレーション方法。
【請求項３】
さらに、
シミュレーショによって計算される対象の前記粒子系は、シミュレーションを行う対象である元粒子系に対して繰り込み処理を行って粒子数を減らしたものであり、
前記元粒子系の粒子に定義される物理量と、くりこまれた前記粒子系の粒子に定義される物理量とは、前記元粒子系とくりこまれた前記粒子系とでレイノルズ数が同一になるように対応付けられている請求項１または２に記載のシミュレーション方法。
【請求項４】
請求項１乃至３のいずれか１項に記載されたシミュレーション方法の手順をコンピュータに実行させるためのプログラム。
【請求項５】
流れ場を形成する複数の粒子を含む粒子系の粒子の各々の位置及び速度を、ある時間刻み幅で分子動力学法により求める際に、
着目している粒子の速度と、着目している粒子の位置における流れ場の流速との差に基づく運動エネルギ、及び着目している粒子の周囲に存在する他の複数の粒子の速度と、前記流速との差に基づく運動エネルギの平均を算出する手順と、
算出された平均値に基づいて、粒子の各々の位置における温度を算出する手順と、
粒子ごとに、各粒子の位置における温度が前記粒子系の目標温度に近づくように温度制御を行うことにより、前記粒子系の温度を目標温度に近づける手順と
を実行するシミュレーション装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、分子動力学法によるシミュレーション方法、プログラム、及びシミュレーション装置に関する。
【背景技術】
【０００２】
分子動力学を用いたシミュレーションが行われている。分子動力学では、シミュレーション対象の系を構成する粒子の運動方程式が数値的に解かれる。分子動力学を用いたシミュレーションを行う場合、シミュレーション対象の系の温度をある目標温度に設定したい場合がある。
【０００３】
下記の特許文献１に、分子動力学によるシミュレーションにおいて、温度設定を変更した場合に、解析状態が迅速に安定化し、所望の設定温度での状態を正確に得ることができるようにする温度制御方法が開示されている。この方法では、分子動力学法によるシミュレーションの途中段階で求められた各粒子の速度を、現時点の系の温度及び目標温度に基づいて更新する。更新された速度に基づいて、系の温度を更新する。この操作により、系の温度を目標温度に近づける。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特開平１０−３３４０７６号公報
【発明の概要】
【発明が解決しようとする課題】
【０００５】
シミュレーション対象の系に流れ場が存在する場合、従来の温度制御方法では、流れ場の流速も温度制御にともなって変化してしまう。このため、流れ場が存在する系のシミュレーション時に従来の温度制御方法を適用すると、シミュレーション結果はシミュレーション対象の系の各粒子の挙動を反映したものにならない。従って、従来の温度制御方法は、流れ場が存在する系のシミュレーションに適用することができない。
【０００６】
本発明の目的は、流れ場が存在する系の分子動力学法によるシミュレーションにおいて、温度制御を行うことが可能なシミュレーション方法を提供することである。
【課題を解決するための手段】
【０００７】
本発明の一観点によると、
流れ場を形成する複数の粒子を含む粒子系を分子動力学法によりシミュレーションする方法であって、
着目している粒子の速度と、着目している粒子の位置における流れ場の流速との差に基づく運動エネルギ、及び着目している粒子の周囲に存在する他の複数の粒子の速度と、前記流速との差に基づく運動エネルギの平均を算出し、
算出された平均値に基づいて、粒子の各々の位置における温度を算出し、
粒子ごとに、各粒子の位置における算出された温度と、前記粒子系の目標温度とに基づいて、温度制御を行うことにより、前記粒子系の温度を目標温度に近づけるシミュレーション方法が提供される。
【０００８】
本発明の他の観点によると、上記シミュレーション方法をコンピュータに実行させるプログラム、及び上記シミュレーション方法を実行するシミュレーション装置が提供される。
【発明の効果】
【０００９】
粒子の速度と、粒子の位置における流れ場の流速との差に基づいて温度制御を行うことにより、流れ場の流速が温度制御の影響を受けにくくなる。これにより、シミュレーション対象の粒子系の流れ場の挙動がシミュレーション結果に反映されることとなる。
【図面の簡単な説明】
【００１０】
【図１】図１は、実施例によるシミュレーション方法のフローチャートである。
【図２】図２Ａは、シミュレーション対象の粒子系の各粒子の位置及び速度を表した模式図であり、図２Ｂは、シミュレーション対象の粒子系の各粒子の位置、及び速度と流速との差を表した模式図である。
【図３】図３は、シミュレーションモデルを示す図である。
【図４】図４Ａ及び図４Ｂは、温度制御を行わない方法を用いてシミュレーションを行った結果を示す図である。
【図５】図５Ａ及び図５Ｂは、流速を考慮しないで温度制御を行う方法を用いてシミュレーションを行った結果を示す図である。
【図６】図６Ａ及び図６Ｂは、実施例による方法でシミュレーションを行った結果を示す図である。
【図７】図７は、実施例によるシミュレーション装置のブロック図である。
【発明を実施するための形態】
【００１１】
図１、図２Ａ〜図２Ｂを参照して、実施例によるシミュレーション方法について説明する。実施例によるシミュレーションにおいて分子動力学法が適用される。分子動力学法では、シミュレーションの対象となる系を複数の粒子で表し、各粒子の運動方程式を数値的に解くことにより粒子の挙動、例えば位置の変化及び速度が求まる。実施例においては、シミュレーション対象である粒子系の複数の粒子によって流れ場が形成される。
【００１２】
図１に、実施例によるシミュレーション方法のフローチャートを示す。以下に説明するシミュレーションは、コンピュータがシミュレーションプログラムを実行することにより実現される。
【００１３】
ステップ１０において、シミュレーション対象の粒子系を定義する物理量、各粒子の初期状態、拘束条件、及び温度制御の実行間隔を設定する。例えば、オペレータの操作によって初期状態、拘束条件、及び温度制御の実行間隔が入力される。粒子系を定義する物理量には、例えば粒子の質量、粒子のポテンシャル等が含まれる。粒子の初期状態には、例えば粒子の位置、粒子の速度、粒子系の温度等が含まれる。粒子系の温度は、粒子系に含まれる複数の粒子の速度の分散に反映される。拘束条件には、境界条件、流れ場の流速、粒子系の目標温度等が含まれる。流れ場の流速は、粒子系に含まれる粒子の速度の平均値に反映される。分子動力学法を用いた動解析の複数のタイムステップのうち、温度制御の実行間隔の設定値で示されたタイムステップ数ごとに温度制御が実行される。
【００１４】
次に、ステップ１１において、現時点で実行するタイムステップが、温度制御を実行するタイムステップに該当するか否かを判定する。例えば、温度制御の実行間隔の設定値で示された回数のタイムステップごとに、温度制御が実行される。
【００１５】
現時点で実行するタイムステップが温度制御を実行するタイムステップに該当しない場合、ステップ１２において、通常の１タイムステップの動解析を実行する。現時点で実行するタイムステップが温度制御を実行するタイムステップに該当する場合、ステップ１３において、１タイムステップの動解析及び温度制御を実行する。
【００１６】
ステップ１２及びステップ１３の後、ステップ１４においてシミュレーションを終了するか否かを判定する。シミュレーションを終了しない場合、ステップ１１に戻ってステップ１１からステップ１４までの処理を繰り返す。シミュレーションを終了する場合には、ステップ１５においてシミュレーション結果を出力する。
【００１７】
次に、ステップ１３の処理について詳細に説明する。
まず、分子動力学法を用いて、１タイムステップの動解析を実行する。１タイムステップの時間刻み幅は、予め設定されている。この動解析により、シミュレーション対象の粒子系の各粒子の位置及び速度が更新される。
【００１８】
１タイムステップの動解析を実行した後、すべての粒子について各粒子の位置における流れ場の流速を算出する。流れ場の流速は、着目する粒子、及び着目する粒子の周囲に存在する他の複数の粒子の速度の平均値で表すことができる。流速を算出する基礎となる複数の粒子として、例えば着目する粒子を含む一部の空間（以下、流速算出単位空間という。）の内側に存在する粒子を採用することができる。流速算出単位空間として、例えば、着目する粒子を中心として、平衡状態のときの粒子間の平均距離の３倍〜４倍程度の半径を持つ球体とすることができる。
【００１９】
次に、図２Ａを参照して、流れ場の流速を算出する具体例について説明する。
図２Ａに、シミュレーション対象の粒子系の各粒子の位置及び速度を表した模式図を示す。図２Ａにおいて、各粒子を丸記号で表し、各粒子の速度ベクトルを矢印で表す。各粒子の速度ベクトルは、流れ場の下流側の方（図２Ａにおいて右方）を向く傾向を示す。
【００２０】
着目する粒子Ｐｉを含む流速算出単位空間３０を定義することができる。着目する粒子Ｐｉの速度をｖｉと表し、流速算出単位空間３０内の他の粒子Ｐｊの速度をｖｊと表す。流速算出単位空間３０内の着目する粒子Ｐｉ以外の粒子の個数をＮＮと表す。着目する粒子Ｐｉの位置における流れ場の流速Ｖｉは、以下の式で表される。
【数１】
0006675785.tif 000002
【００２１】
流速Ｖｉを算出した後、各粒子Ｐｉ、Ｐｊの速度ｖｉ、ｖｊと、着目する粒子Ｐｉの位置における流れ場の流速Ｖｉとの差に基づいて、着目する粒子Ｐｉの位置における温度を算出する。
【００２２】
図２Ｂに、シミュレーション対象の粒子系の各粒子Ｐｉ、Ｐｊの位置、及び速度ｖｉ、ｖｊと流速Ｖｉとの差を表した模式図を示す。各粒子に付された矢印が、その粒子の速度ｖｉ、ｖｊと流速Ｖｉとの差を表している。粒子の速度ｖｉ、ｖｊと流速Ｖｉとの差を流速算出単位空間３０内の粒子について平均した速度は０になる。
【００２３】
着目する粒子Ｐｉの位置における温度をＴｉ、ボルツマン定数をｋＢ、粒子Ｐｉ、Ｐｊの質量をそれぞれｍｉ、ｍｊと表すと、以下の式（２）が成り立つ。式（２）を用いて、着目する粒子Ｐｉの位置における温度Ｔｉを算出することができる。
【数２】
0006675785.tif 000003
【００２４】
式（２）の右辺は、着目している粒子Ｐｉの速度ｖｉと流速Ｖｉとの差に基づく運動エネルギ、及び着目している粒子Ｐｉの周囲に存在する他の複数の粒子Ｐｊの速度ｖｊと流速Ｖｉとの差に基づく運動エネルギの平均値を表している。この平均値に基づいて、着目する粒子Ｐｉの位置における温度Ｔｉが算出される。
【００２５】
温度Ｔｉを算出した後、各粒子の位置における温度が目標温度に近づくように、各粒子の速度を調整する。速度を調整する一手法として、速度スケーリング法を適用することができる。速度スケーリング法においては、着目する粒子Ｐｉの速度調整後の速度をｖｉｍと表すと、速度ｖｉｍは、以下の式で表される。
【数３】
0006675785.tif 000004
ここで、ａはスケーリングの強さを決めるパラメータであり、１以下の正の実数である。ｓｉは、スケーリングファクタと呼ばれる。Ｔｃは目標温度である。
【００２６】
粒子の速度を調整することにより、１タイムステップの動解析及び温度制御の処理が終了する。
【００２７】
次に、分子動力学法によるシミュレーションにおいて、流れ場の流速を考慮しないで温度制御を適用する方法と比較して、上記実施例の効果について説明する。流れ場の流速を考慮しないで温度制御を行う場合には、例えば、着目する粒子Ｐｉの位置における温度Ｔｉは上述の式（２）に代えて、以下の式（４）を用いて算出される。
【数４】
0006675785.tif 000005
【００２８】
さらに、温度制御による速度調整は、式（３）に代えて、以下の式（５）が適用される。
【数５】
0006675785.tif 000006
【００２９】
上述の式（４）の速度ｖｉ、ｖｊには、図２Ａに示したように流れ場の流速Ｖｉが含まれている。流速Ｖｉを含んだ速度ｖｉ、ｖｊに基づいて温度制御が行われるため、温度制御によって流れ場の流速Ｖｉが乱されてしまう。上記実施例による方法では、式（３）に示したように、粒子の速度ｖｉ、ｖｊと、粒子Ｐｉの位置における流れ場の流速Ｖｉとの差に基づいて温度制御が行われ、温度制御実施後の速度に流速Ｖｉが加算される。このため、温度制御による流速Ｖｉの乱れを抑制することができる。
【００３０】
上記実施例の効果を確認するために、同一のシミュレーションモデルを用いて、温度制御を行わない方法、流速を考慮しないで温度制御を行う方法、及び実施例による方法の３通りの方法でシミュレーションを行った。次に、これらのシミュレーション結果について説明する。
【００３１】
図３にシミュレーションモデルを示す。シミュレーションは２次元空間で行った。２次元の正方形の領域内に円形のシリンダが配置されている。この空間に流れる流体を構成する粒子の挙動をシミュレーションにより求めた。ｘｙ直交座標系を定義し、正方形の領域の１つの頂点を原点とし、正方形の各辺を、ｘ軸またはｙ軸に平行にした。正方形の一辺の長さＬ１を３９２２×１０−１０ｍとし、シリンダの直径Ｄを５００×１０−１０ｍとした。シリンダは、ｙ方向に関して正方形の領域の中心に配置されている。正方形の領域のｘ方向の負側の辺（ｘ＝０）からシリンダの中心までの距離Ｌ２をＬ１／４＝９８０．５×１０−１０ｍとした。
【００３２】
流体を構成する粒子の質量を６．６３×１０−２６ｋｇとした。正方形の領域内に配置される粒子の個数は約２６万個である。粒子間の相互作用ポテンシャルとして、レナードジョーンズポテンシャルの斥力の項のみを用いた。具体的には、相互作用ポテンシャルφ（ｒ）は、以下の式で定義される。
【数６】
0006675785.tif 000007
ここで、ε＝１．６５×１０−２１Ｊ、σ＝３．４１×１０−１０ｍである。この粒子系の流体の粘度は３．６４×１０−４Ｐａ・ｓであり、密度は１．１８×１０３ｋｇ／ｍ３である。
【００３３】
さらに、シリンダは、円周上に相互にバネで接続された複数の粒子を配置することにより表現した。シリンダ表面に配置した粒子に対しては、ランジュバン法による温度制御を行った。シリンダの温度は８８．８５Ｋとした。
【００３４】
温度制御を行う場合の目標温度Ｔｃとして、温度制御を行わない方法でシミュレーションしたときの定常状態到達時の平均温度である２３５．１５Ｋを用いた。
【００３５】
次に、境界条件について説明する。右端（ｘ＝Ｌ１）と左端（ｘ＝０）とで特殊な周期境界条件が課されている。具体的には、粒子が右端の境界を通過するとき、ｘ方向の流速６００ｍ／ｓを中心とし、温度を８８．８５Ｋとするボルツマン分布に基づいて、ランダムに新しい速度が割り当てられる。この粒子が、左端の境界から正方形の領域内に流入する。下端（ｙ＝０）と上端（ｙ＝Ｌ１）の境界は、ｘ方向の流速が６００ｍ／ｓ、温度が８８．８５Ｋの熱浴として働く。下端及び上端の境界に到達した粒子は、ｘ方向の流速６００ｍ／ｓだけシフトしたハーフマクスウェル分布に基づいて、ランダムに新しい速度が割り当てられる。
【００３６】
図４Ａ、図４Ｂ、図５Ａ、図５Ｂ、図６Ａ、及び図６Ｂにシミュレーショ結果を示す。図４Ａ及び図４Ｂは、温度制御を行わない方法を用いてシミュレーションを行った結果を示し、図５Ａ及び図５Ｂは、流速を考慮しないで温度制御を行う方法を用いてシミュレーションを行った結果を示し、図６Ａ及び図６Ｂは、実施例による方法でシミュレーションを行った結果を示す。
【００３７】
図４Ａ、図５Ａ、及び図６Ａは、流体に生じた渦を可視化した図である。渦を可視化する方法は、例えば特許第５８８８９２１号公報に開示されている。図４Ｂ、図５Ｂ、及び図６Ｂは、温度分布を濃淡で表す図である。色の相対的に濃い領域が相対的に高温であることを表している。図４Ａ〜図６Ｂまでの各々において、左から右に並ぶ複数の図は、時間の経過に対応している。
【００３８】
図４Ａに示されているように、温度制御を行わない場合には、シリンダの後方にカルマン渦が発生していることを確認できる。また、図４Ｂに示されているように、シリンダ近傍を通過する流体の温度が相対的に大きく上昇しおり、顕著な温度分布が発生していることがわかる。定常状態に達した時点における流体の平均温度は２３５．１５Ｋであった。このように、顕著な温度分布が発生したのは、シミュレーションにおいて温度制御を行っていないためである。
【００３９】
流速を考慮しないで温度制御を行った場合には、図５Ｂに示されているように、流体の温度は、図４Ｂの場合に比べて均一になっている。これは、温度制御が行われたことの効果である。ところが、図５Ａに示されているように、カルマン渦が確認できない。これでは、流体の挙動が適切にシミュレーションされたとはいえない。これは、温度制御時に流速を考慮していないため、温度制御が流速自体に影響を与えてしまったためである。
【００４０】
図６Ａに示されているように、実施例による方法で温度制御を行った場合には、カルマン渦が確認される。さらに、図６Ｂに示されているように、流体の温度がほぼ均一になっている。実施例による方法を採用することにより、流れ場を大きく乱すことなく、かつ適切な温度制御を行うことが可能であることが確認された。
【００４１】
次に、上記実施例の変形例について説明する。上記実施例では、式（３）を用いて粒子の速度ｖｉを調整することにより温度制御を行った。その他に、速度に比例する摩擦力を与えて温度制御を行ってもよいし、Ｎｏｓｅ−Ｈｏｏｖｅｒ法を用いて温度制御を行ってもよい。
【００４２】
速度に比例する摩擦力を与えて温度制御を行う場合には、ステップ１３（図１）において動解析を行う際に、通常の運動方程式に代えて、以下の式（７）を用いることができる。
【数７】
0006675785.tif 000008
ここで、ｒｉは着目する粒子Ｐｉの位置、ｆｉは粒子Ｐｉに作用する力を表す。γはスケーリングの強さを決めるパラメータである。
【００４３】
Ｎｏｓｅ−Ｈｏｏｖｅｒ法を用いて温度制御を行う場合には、通常の運動方程式に代えて、以下の式（８）を用いることができる。
【数８】
0006675785.tif 000009
ここで、ｒｉは粒子Ｐｉの位置、Φは粒子間ポテンシャル、ｇは空間の自由度を表す。Ｑはスケーリングの強さを決めるパラメータである。
【００４４】
分子動力学法によるシミュレーション方法において、計算対象の粒子数を減らすために、種々の方法が提案されている。例えば、流路の代表長さを短くする方法、繰り込みの手法を用いて粒子数を減らす方法（例えば、特許第５２４１４６８号公報）等により粒子数を減らすことができる。
【００４５】
シミュレーション対象である元粒子系の粒子数を減らすとき、シミュレーション対象である元粒子系と、粒子数を減らした粒子系とで、レイノルズ数が同一になるように粒子に与えられる物理量を対応付けることが好ましい。レイノルズ数を同一にすることにより、元粒子系と粒子数を減らした粒子系との間で、流れの相似性を確保することができる。
【００４６】
レイノルズ数Ｒｅは、以下の式（９）で与えられる。
【数９】
0006675785.tif 000010
ここで、μは流体の粘度を表し、ρは流体の密度を表し、ｖは粒子の速度を表し、Ｌは流路の代表長さを表す。例えば、粒子数を減らすために流路を小さくして代表長さＬを１／ｘ倍にした場合、粒子数を減らした粒子系の粒子の速度ｖをｘ倍にすればよい。
【００４７】
繰り込みの手法を用いて粒子数を減らす場合には、例えば、繰り込みによって粒子系の密度ρ、代表長さＬが変化せず、粘度μがλ倍になる場合がある。ここで、λは繰り込み因子であり、繰り込み回数ｎを用いてλ＝２ｎと表される。この場合には、レイノルズ数Ｒｅを不変とするために、速度ｖをλ倍にすればよい。
【００４８】
図７に、実施例によるシミュレーション装置のブロック図を示す。このシミュレーション装置として、シミュレーションプログラムを実行するコンピュータを用いることができる。このコンピュータは、中央処理装置（ＣＰＵ）４０、メインメモリ４１、補助記憶装置４２、入力装置４３、及び出力装置４４を含む。補助記憶装置４２は、シミュレーションプログラムが格納されている記録媒体を含む。この記録媒体は、補助記憶装置に内蔵されたものでもよく、補助記憶装置に対して着脱可能なリムーバブル媒体でもよい。このシミュレーションプログラムがメインメモリ４１にロードされ、ＣＰＵ４０によって実行される。
【００４９】
入力装置４３からシミュレーションに必要な情報が入力される。例えば、ステップ１０（図１）で設定される種々の情報が入力される。シミュレーション結果が出力装置４４に出力される。例えば、ステップ１５（図１）において、図６Ａに示した流体に生じた渦を可視化した画像、図６Ｂに示した流体の温度分布を示す画像等が出力装置４４に表示される。
【００５０】
図７に示したシミュレーション装置により、上記実施例によるシミュレーションを行うことができる。
【００５１】
上述の実施例及び変形例は例示であり、異なる実施例及び変形例で示した構成の部分的な置換または組み合わせが可能であることは言うまでもない。複数の実施例及び変形例の同様の構成による同様の作用効果については実施例及び変形例ごとには逐次言及しない。さらに、本発明は上述の実施例及び変形例に制限されるものではない。例えば、種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。
【符号の説明】
【００５２】
１０〜１５  ステップ
３０  流速算出単位空間
４０  中央処理装置（ＣＰＵ）
４１  メインメモリ
４２  補助記憶装置
４３  入力装置
４４  出力装置
【図１】
0006675785.tif 000011
【図２】
0006675785.tif 000012
【図３】
0006675785.tif 000013
【図４】
0006675785.tif 000014
【図５】
0006675785.tif 000015
【図６】
0006675785.tif 000016
【図７】
0006675785.tif 000017