(57)【要約】【課題】高レイノルズ数の条件下での粒子の挙動のシミュレーションに適したシミュレーション方法を提供する。【解決手段】シミュレーション対象の粒子系の粒子間の相互作用ポテンシャルを、フィッティングパラメータを含む関数で定義して、定義された関数に基づいて粒子系の挙動を分子動力学法によりシミュレーションする。このとき、粒子系の粘度が低下し、かつ粒子系における音速が不変となる条件を満たすようにフィッティングパラメータを変換する。変換後のフィッティングパラメータを持つ関数を用いて、分子動力学法によるシミュレーションを実行する。【選択図】図１
2018049583.tif 000002
【特許請求の範囲】
【請求項１】
シミュレーション対象の粒子系の粒子間の相互作用ポテンシャルを、フィッティングパラメータを含む関数で定義して、定義された前記関数に基づいて前記粒子系の挙動を分子動力学法によりシミュレーションする方法であって、
前記粒子系の粘度が低下し、かつ前記粒子系における音速が不変となる条件を満たすように前記フィッティングパラメータを変換し、
変換後の前記フィッティングパラメータを持つ前記関数を用いて、分子動力学法によるシミュレーションを実行するシミュレーション方法。
【請求項２】
前記粒子系の粒子の各々に作用するケージポテンシャルを、着目する粒子の周囲の粒子に起因する前記関数の和で定義したとき、前記ケージポテンシャルの極値における二次導関数の値が、前記フィッティングパラメータの変換の前後で不変である請求項１に記載のシミュレーション方法。
【請求項３】
粒子間の距離ｒに依存する前記関数φ（ｒ）は、エネルギの次元を持つフィッティングパラメータε、長さの次元を持つフィッティングパラメータσ、平衡時の粒子間の距離ｒ０、及び無次元関数ｆを用いて、
【数１】
2018049583.tif 000024
と表される請求項１または２に記載のシミュレーション方法。
【請求項４】
前記関数は、隣り合う粒子が斥力を及ぼし合う斥力領域と、引力を及ぼし合う引力領域とに分けられ、斥力領域においては傾きの絶対値が大きくなり、引力領域においては傾きの絶対値が小さくなるように、斥力領域と引力領域とで異なる変換則に基づいて前記フィッティングパラメータを変換する請求項１乃至３のいずれか１項に記載のシミュレーション方法。
【請求項５】
前記関数は、隣り合う粒子が斥力を及ぼし合う斥力領域と、引力を及ぼし合う引力領域とに分けられ、引力領域においてはポテンシャルの深さが浅くなり、かつポテンシャルの幅が狭くなるように前記フィッティングパラメータを変換し、斥力領域においてはポテンシャルの形状が変化せず、かつ引力領域のポテンシャルに連続するようにポテンシャルの深さを浅くする請求項１乃至３のいずれか１項に記載のシミュレーション方法。
【請求項６】
請求項１乃至５のいずれか１項に記載のシミュレーション方法をコンピュータで実行するためのコンピュータプログラム。
【請求項７】
シミュレーション対象の粒子系の粒子間の相互作用ポテンシャルを、フィッティングパラメータを含む関数で定義して、定義された前記関数に基づいて前記粒子系の挙動を分子動力学法によりシミュレーションする処理装置を含み、
前記処理装置は、前記粒子系の粘度が小さくなり、かつ前記粒子系における音速が不変となる条件を満たすように、前記関数の前記フィッティングパラメータを、変換係数に基づいて変換する変換則を記憶しており、
入力された前記変換係数、及び前記変換則に基づいて前記フィッティングパラメータを変換し、
変換後の前記フィッティングパラメータを持つ前記関数を用いて、分子動力学法によるシミュレーションを実行するシミュレーション装置。
【請求項８】
さらに、前記フィッティングパラメータ、及び前記変換係数を入力する入力装置を有する請求項７に記載のシミュレーション装置。
【請求項９】
前記粒子系の粒子の各々に作用するケージポテンシャルを、着目する粒子の周囲の粒子に起因する前記関数の和で定義したとき、前記ケージポテンシャルの極値における二次導関数の値が、前記フィッティングパラメータの変換の前後で不変である請求項７または８に記載のシミュレーション装置。
【請求項１０】
粒子間の距離ｒに依存する前記関数φ（ｒ）は、エネルギの次元を持つフィッティングパラメータε、長さの次元を持つフィッティングパラメータσ、平衡時における粒子間の距離を規定する定数ｒ０、及び無次元関数ｆを用いて、
【数２】
2018049583.tif 000025
と表される請求項７乃至９のいずれか１項に記載のシミュレーション装置。
【請求項１１】
前記関数は、隣り合う粒子が斥力を及ぼし合う斥力領域と、引力を及ぼし合う引力領域とに分けられ、
前記変換則は、前記関数の斥力領域と引力領域とで異なっており、斥力領域においては前記相互作用ポテンシャルの傾きの絶対値が大きくなり、引力領域においては傾きの絶対値が小さくなるように定義されている請求項７乃至１０のいずれか１項に記載のシミュレーション装置。
【請求項１２】
前記関数は、隣り合う粒子が斥力を及ぼし合う斥力領域と、引力を及ぼし合う引力領域とに分けられ、
前記変換則は、引力領域においてはポテンシャルの深さが浅くなり、かつポテンシャルの幅が狭くなるように前記フィッティングパラメータを変換し、斥力領域においてはポテンシャルの形状が変化せず、かつ引力領域のポテンシャルに連続するように定義されている請求項７乃至１１のいずれか１項に記載のシミュレーション装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、分子動力学法を用いて粒子系の挙動をシミュレーションするシミュレーション方法、コンピュータプログラム、及びシミュレーション装置に関する。
【背景技術】
【０００２】
分子動力学法を用いたコンピュータシミュレーションが行われている。分子動力学法では、シミュレーションの対象となる系を構成する粒子の運動方程式が数値的に解かれる。シミュレーション対象の系に含まれる粒子数が増加すれば、必要な計算量が増大する。現行のコンピュータの演算能力でシミュレーション可能な系の粒子数は、通常、数十万個程度である。
【０００３】
下記の特許文献１に、くりこみ変換の手法を用いたシミュレーション方法が開示されている。くりこみ変換されて粒子数が減少した粒子系について分子動力学法による計算を行うことにより、シミュレーションに必要な計算量を減らすことができる。
【０００４】
特許文献１に開示された方法でくりこみ変換された粒子系の寸法は、元の粒子系の寸法より小さい。この手法を流体のシミュレーションに適用する場合、くりこみ変換によって粒子系の寸法が小さくなると、流路の代表的な長さも短くなる。流路の代表的な長さの低下は、レイノルズ数の低下に繋がる。
【先行技術文献】
【特許文献】
【０００５】
【特許文献１】特許第５２４１４６８号公報
【発明の概要】
【発明が解決しようとする課題】
【０００６】
粒子系に対してくりこみ変換処理を行うと、レイノルズ数を低下させることになる。従来のくりこみ群分子動力学法では、高レイノルズ数の条件下での粒子の挙動をシミュレーションすることが困難である。
【０００７】
本発明の目的は、高レイノルズ数の条件下での粒子の挙動のシミュレーションに適したシミュレーション方法、コンピュータプログラム、及びシミュレーション装置を提供することである。
【課題を解決するための手段】
【０００８】
本発明の一観点によると、
シミュレーション対象の粒子系の粒子間の相互作用ポテンシャルを、フィッティングパラメータを含む関数で定義して、定義された前記関数に基づいて前記粒子系の挙動を分子動力学法によりシミュレーションする方法であって、
前記粒子系の粘度が低下し、かつ前記粒子系における音速が不変となる条件を満たすように前記フィッティングパラメータを変換し、
変換後の前記フィッティングパラメータを持つ前記関数を用いて、分子動力学法によるシミュレーションを実行するシミュレーション方法が提供される。
【０００９】
本発明の他の観点によると、
上記シミュレーション方法をコンピュータで実行するためのコンピュータプログラムが提供される。
【００１０】
本発明のさらに他の観点によると、
シミュレーション対象の粒子系の粒子間の相互作用ポテンシャルを、フィッティングパラメータを含む関数で定義して、定義された前記関数に基づいて粒子系の挙動を分子動力学法によりシミュレーションする処理装置を含み、
前記処理装置は、前記粒子系の粘度が小さくなり、かつ前記粒子系における音速が不変となる条件を満たすように、前記関数の前記フィッティングパラメータを、変換係数に基づいて変換する変換則を記憶しており、
入力された前記変換係数、及び前記変換則に基づいて前記フィッティングパラメータを変換し、
変換後の前記フィッティングパラメータを持つ前記関数を用いて、分子動力学法によるシミュレーションを実行するシミュレーション装置が提供される。
【発明の効果】
【００１１】
音速が不変となる条件の下で粒子系の粘度を低下させることにより、粒子の速度を低下させることなく、粘度を低下させた条件の下でシミュレーションを行うことが可能になる。レイノルズ数は、粒子の速度に比例し、粘度に反比例する。従って、高レイノルズ数の条件で、分子動力学法を用いて流体の挙動をシミュレーションすることが可能になる。
【図面の簡単な説明】
【００１２】
【図１】図１は、実施例によるシミュレーション方法のフローチャートである。
【図２】図２は、相互作用ポテンシャルを表す関数φ（ｒ）の一例を示すグラフである。
【図３】図３は、ケージポテンシャルの例を示すグラフである。
【図４】図４Ａは、シミュレーションにより求めた動粘度と定数ｂとの関係を示すグラフであり、図４Ｂは、定数ｂを変化させたときの各パラメータの値を示す図表である。
【図５】図５は、他の実施例による相互作用ポテンシャルを表す関数φ（ｒ）の形状を示すグラフである。
【図６】図６は、他の実施例によるシミュレーションにより求めた動粘度と定数ｂとの関係を示すグラフである。
【図７】図７は、実施例によるシミュレーション装置の機能ブロック図である。
【発明を実施するための形態】
【００１３】
分子動力学法によるシミュレーションでは、粒子間の相互作用ポテンシャルとして、レナードジョーンズポテンシャル、モースポテンシャル等が用いられる。これらのポテンシャルは、フィッティングパラメータε、σ、及び無次元関数ｆを用いて以下の関数φ（ｒ）で表すことができる。
【数１】
2018049583.tif 000003
式（１）中のフィッティングパラメータεはエネルギの次元を持ち、フィッティングパラメータσは長さの次元を持つ。ｒは粒子間の距離を表し変数であり、ｒ０は平衡時の粒子間距離を表す定数である。
【００１４】
レイノルズ数Ｒｅは、以下の式（２）で表される。
【数２】
2018049583.tif 000004
ここで、μは流体の粘度を表し、ρは流体の密度を表し、ｖは粒子の速度を表し、Ｌは流路の代表長さを表す。本願発明者の評価実験によると、相互作用ポテンシャルを定義する関数φ（ｒ）のフィッティングパラメータε及びσを調節して粒子系の粘度μを変更すると、粒子系における音速ｃが粘度μに比例して変化することが判明した。流体のシミュレーションは、各粒子の速度ｖが音速ｃを超えない範囲で行うことが好ましい。
【００１５】
高レイノルズ数の条件下、すなわち乱流が発生する条件下でシミュレーションを行うには、式（２）から粘度μを小さくすればよいことが分かる。ところが、粘度μを小さくすると音速ｃもそれに比例して小さくなってしまうため、粒子の速度ｖの取り得る範囲の上限値も小さくなる。粒子の速度ｖの範囲の上限値が小さくなることは、レイノルズ数Ｒｅを高くすることに反する。従って、単にフィッティングパラメータε及びσを調節して粘度を変更する方法では、レイノルズ数Ｒｅの高い条件でシミュレーションを行うことが困難である。
【００１６】
以下に説明する実施例においては、粒子系の粘度が小さくなり、かつ音速は不変となる条件を満たすように、フィッティングパラメータε及びσを変換し、関数φ（ｒ）のフィッティングパラメータε及びσを変換後の値に置き換えた関数を用いて分子動力学法による計算を行う。
【００１７】
次に、図１〜図３を参照して、実施例によるシミュレーション方法について説明する。
【００１８】
図１は、実施例によるシミュレーション方法のフローチャートである。まず、シミュレーション対象の粒子系の粒子間の相互作用ポテンシャルを表す関数φ（ｒ）を定義する（ステップＳ１）。関数φ（ｒ）は、例えば、上述の式（１）のように定義することができる。
【００１９】
図２は、相互作用ポテンシャルを表す関数φ（ｒ）の一例を示すグラフである。横軸は粒子の中心からの距離ｒを表し、縦軸は関数φ（ｒ）の値を表す。関数φ（ｒ）を太い実線で示している。関数φ（ｒ）はｒ＝ｒ０において極小値を示す。距離ｒ０は、平衡時における粒子間距離に相当する。距離ｒがｒ０以下の領域では、粒子間に斥力が作用し、距離ｒがｒ０以上の領域では、粒子間に引力が作用する。距離ｒがｒ０以下の領域を「斥力領域」といい、距離ｒがｒ０以上の領域を「引力領域」ということとする。
【００２０】
関数φ（ｒ）を定義した後、斥力領域と引力領域とで異なる変換則に基づいて、フィッティングパラメータε及びσを変換することにより、斥力領域におけるフィッティングパラメータε１及びσ１、引力領域におけるフィッティングパラメータε２及びσ２を求める（ステップＳ２）。斥力領域における相互作用ポテンシャルを表す関数φ１（ｒ）及び引力領域における相互作用ポテンシャルを表す関数φ２（ｒ）は、以下の式で表される。
【数３】
2018049583.tif 000005
【００２１】
図２に、関数φ（ｒ）とともに、関数φ１（ｒ）及びφ２（ｒ）の一例を破線で示す。関数φ１（ｒ）とφ２（ｒ）とは、ｒ＝ｒ０における値が等しく、かつ傾きは共に０である。また、距離ｒが無限大の極値において、関数φ２（ｒ）の値は、関数φ（ｒ）の値に等しい。斥力領域においては、関数φ１（ｒ）の傾きの絶対値が関数φ（ｒ）の絶対値より大きくなっている。これに対し、引力領域では、関数φ２（ｒ）の傾きの絶対値が関数φ（ｒ）の傾きの絶対値より小さくなっている。
【００２２】
フィッティングパラメータε及びσの変換を行なった後、関数φ１（ｒ）及びφ２（ｒ）で定義される相互作用ポテンシャルを用いて分子動力学法により粒子系の挙動のシミュレーションを行う（ステップＳ３）。
【００２３】
次に、フィッティングパラメータ変換後の関数φ１（ｒ）及び関数φ２（ｒ）で定義される相互作用ポテンシャルを持つ粒子系の粘度について説明する。粒子系の粘度（動粘度）は、周知のＳＬＬＯＤ法を用いて求めることができる。ＳＬＬＯＤ法については、Evans, D.J., Hoover, W.G., Failor, B.H., Moran, B., and Ladd, A.J.C., Phys. Rev. A28 (1983), 1016に説明されている。
【００２４】
図２に示した関数φ１（ｒ）及びφ２（ｒ）で定義される相互作用ポテンシャルについてＳＬＬＯＤ法を用いて粘度を求めたところ、求められた粘度は、関数φ（ｒ）で定義される相互作用ポテンシャルについて求めた粘度より低くなっていることが確認された。
【００２５】
次に、フィッティングパラメータ変換後の関数φ１（ｒ）及び関数φ２（ｒ）で定義される相互作用ポテンシャルを持つ粒子系の音速について説明する。
【００２６】
隣り合う粒子をばねで連結した粒子系のヤング率Ｅは、ばね定数ｋに比例することが知られている。また、音速ｃは、ヤング率Ｅの平方根に比例することが知られている。すなわち、音速ｃは、ばね定数ｋの平方根に比例する。粒子系のばね定数ｋは、各粒子に作用するケージポテンシャルの極値における二次導関数の値で与えられる。このため、音速ｃは、ケージポテンシャルの極値における二次導関数の値の平方根に比例する。従って、ケージポテンシャルの極値における二次導関数の値が、フィッティングパラメータの変換の前後において不変であれば、音速ｃも不変であるといえる。
【００２７】
次に、図３を参照してケージポテンシャルについて説明する。
【００２８】
図３は、ケージポテンシャルの例を示すグラフである。図３の横軸は着目する粒子からの距離ｒを表し、縦軸はポテンシャルの大きさを表す。曲線ａ１は、ｒ＝−ｒ０の位置の粒子に起因する相互作用ポテンシャルを表し、破線ａ２は、ｒ＝ｒ０の位置の粒子に起因する相互作用ポテンシャルを表す。ｒ＝０の位置の粒子に作用するケージポテンシャルφｃは、この曲線ａ１と破線ａ２で表される２つの相互作用ポテンシャルの和で表される。曲線ａ１は、図２に示した相互作用ポテンシャルを左方向にｒ０だけ平行移動させることにより得られる。破線ａ２は、ｒ＝０に関して曲線ａ１と線対称である。
【００２９】
次に、粒子系の相互作用ポテンシャルが、フィッティングパラメータε及びσの変換を行う前の式（１）で表された関数φ（ｒ）で定義される場合の音速ｃについて説明する。ｒ＝−ｒ０の位置の粒子による相互作用ポテンシャルを表す関数φｎ（ｒ）、及びｒ＝ｒ０の位置の粒子による相互作用ポテンシャルを表す関数φｐ（ｒ）は、それぞれ以下の式で表される。
【数４】
2018049583.tif 000006
【００３０】
ケージポテンシャルφｃ（ｒ）は、以下の式で表される。
【数５】
2018049583.tif 000007
【００３１】
ケージポテンシャルφｃ（ｒ）の二次導関数は、以下の式で表される。
【数６】
2018049583.tif 000008
【００３２】
音速ｃは、ケージポテンシャルの極値における二次導関数の値の平方根に比例するため、以下の式が導出される。
【数７】
2018049583.tif 000009
【００３３】
次に、粒子系の相互作用ポテンシャルが、フィッティングパラメータε及びσを変換した後の式（３）で表された関数φ１（ｒ）及びφ２（ｒ）で定義される場合の音速について説明する。ｒ≦０の領域のケージポテンシャルφｃｎ（ｒ）及びｒ≧０の領域のケージポテンシャルφｃｐ（ｒ）は、それぞれ以下の式で表される。
【数８】
2018049583.tif 000010
ケージポテンシャルφｃｎ（ｒ）とφｃｐ（ｒ）とは、原点に関して相互に対称であるため、以下、一方のケージポテンシャルφｃｎ（ｒ）について説明する。
【００３４】
ケージポテンシャルφｃｎ（ｒ）の二次導関数は以下の式で表される。
【数９】
2018049583.tif 000011
【００３５】
音速ｃは、ケージポテンシャルの極値における二次導関数の値の平方根に比例するため、以下の式が導出される。
【数１０】
2018049583.tif 000012
【００３６】
図３では、着目する粒子の両側の粒子に起因する相互作用ポテンシャルを表す関数の和によりケージポテンシャルを定義したが、二次元空間または三次元空間においては、着目する粒子の周囲に存在する粒子に起因する相互作用ポテンシャルを表す関数の和によってケージポテンシャルが定義される。
【００３７】
式（７）と式（１０）とを比較すると、フィッティングパラメータの変換の前後で音速が不変となるためには、以下の条件を満たせばよいことがわかる。
【数１１】
2018049583.tif 000013
【００３８】
次に、上述の式（１１）を満たすようにフィッティングパラメータの変換を行った具体例について説明する。式（１１）を満足する一例として、定数（変換係数）ｂを用いた下記の変換則を適用する。
【数１２】
2018049583.tif 000014
κは、ケージポテンシャルを構成する際に、引力領域による音速の低下を補正する。定数ｂを無限大にした極限が、ファンデルワールス近似に対応する。
【００３９】
式（４）〜式（１１）で説明したように、式（１２）の変換則を用いてフィッティングパラメータε及びσを変換した後のケージポテンシャルの極値における二次導関数の値は、変換前のケージポテンシャルの極値における二次導関数の値と等しく、定数ｂに依存しない。このため、フィッティングパラメータε及びσの変換前後において、粒子系の音速は不変である。
【００４０】
次に、図４を参照して、式（１２）の変換則を適用してフィッティングパラメータε及びσを変換して得られた相互作用ポテンシャルを持つ粒子系の粘度について説明する。
【００４１】
粒子系の相互作用ポテンシャルとしてレナードジョーンズ型ポテンシャルを適用し、式（１２）の定数ｂを変化させて、ＳＬＬＯＤ法を用いて粒子系の動粘度をシミュレーションにより求めた。
【００４２】
図４Ａは、シミュレーションにより求めた動粘度と定数ｂとの関係を示すグラフである。横軸は定数ｂを表し、縦軸は動粘度を単位「ｍ２／ｓ」で表す。定数ｂ＝１のとき、すなわちフィッティングパラメータε及びσの変換前において、シミュレーションの条件を下記の通りとした。
フィッティングパラメータε＝４０４．５２Ｋ
フィッティングパラメータσ＝２．６４×１０−１０ｍ
平衡時の粒子間の距離ｒ０＝２．６４×１０−１０ｍ
粒子の質量ｍ＝１８．０２Ｄａ
粒子系の温度Ｔ＝３００．０Ｋ
粒子系の圧力Ｐ＝２０．０ａｔｍ（約２．０ＭＰａ）
密度ρ＝１１８５．４ｋｇ／ｍ３
【００４３】
図４Ｂに、定数ｂを変化させたときの各パラメータの値を示す。温度の次元を持つフィッティングパラメータεを調整したときに無次元温度ｋＢＴ／εが一定になるように温度Ｔを調整している。無次元温度が一定であれば、シミュレーション結果は等しくなる。密度ρの変化は、極値近傍のポテンシャルの非対称性によって粒子間距離がずれることにより生じている。定数ｂを変化させたことによる密度ρの変化は、動粘度の変化に比べて極わずかである。
【００４４】
図４Ａに示したシミュレーション結果からわかるように、定数ｂを大きくすると動粘度が低下する。例えば、定数ｂを３２〜６４程度にすると、動粘度がフィッティングパラメータ変換前の値に比べて１／１０程度まで低下している。
【００４５】
図４Ａ及び図４Ｂに示したシミュレーション結果から、式（１１）を満たすようにフィッティングパラメータε及びσを変換することにより、音速を一定に維持したまま、動粘度を低下させることができることが確認された。
【００４６】
次に、上記実施例の優れた効果について説明する。
上述の式（３）及び式（１１）に基づいて、相互作用ポテンシャルを定義する関数を決定すると、粒子系の音速を一定に維持したまま、動粘度を低下させることができる。これにより、高レイノルズ数の条件下で流体の挙動を分子動力学法によりシミュレーションすることが可能になる。
【００４７】
次に、図５及び図６を参照して、他の実施例によるシミュレーション方法について説明する。以下、上記実施例と共通の構成については説明を省略する。
【００４８】
本実施例においては、上記実施例の式（３）に代えて、斥力領域における相互作用ポテンシャルを表す関数φ１（ｒ）、及び引力領域における相互作用ポテンシャルを表す関数φ２（ｒ）を、以下のように定義する。
【数１３】
2018049583.tif 000015
ここで、Ｃは、ｒ＝ｒ０で関数φ（ｒ）が連続するように、すなわちφ１（ｒ０）＝φ２（ｒ０）が成立するように調整するための調整定数である。本実施例では、図１のステップＳ２において、式（１３）に基づいて関数φ１（ｒ）及びφ２（ｒ）を決定する。
【００４９】
式（１３）の斥力領域における相互作用ポテンシャルを表す関数φ１（ｒ）の定義式の第２項は定数であるため、音速ｃを規定する式（７）の二次導関数には影響を与えない。このため、本実施例においても、音速不変の条件式（１１）を適用することができる。本実施例では、条件式（１１）を満たす変換則として、下記の変換則を適用する。
【数１４】
2018049583.tif 000016
【００５０】
以下、相互作用ポテンシャルを表す関数φ（ｒ）として、以下に示すシフティドレナードジョーンズポテンシャルを用いた例について説明する。
【数１５】
2018049583.tif 000017
【００５１】
式（１５）で定義したシフティドレナードジョーンズポテンシャルを表す関数φ（ｒ）を式（１）に当てはめると、関数ｆ（ｘ）は以下の式で表される。
【数１６】
2018049583.tif 000018
【００５２】
一般的なレナードジョーンズポテンシャルにおいては、フィッティングパラメータσ（ポテンシャルの幅）の変更に伴いｒ０（ポテンシャルの鞍点）が変わってしまう。一般的なレナードジョーンズポテンシャルでは、以下の式が成立する。
【数１７】
2018049583.tif 000019
このため、フィッティングパラメータσが決まると、ｒ０も自動的に決まってしまう。これに対し、式（１５）のシフティドレナードジョーンズポテンシャルにおいては、フィッティングパラメータσ（ポテンシャルの幅）とｒ０（ポテンシャルの鞍点）とを独立に設定することが可能になる。
【００５３】
相互作用ポテンシャルを表す関数φを式（１５）で定義し、式（１４）のフィッティングパラメータの変換則を適用した場合、ｒ＝ｒ０でφ１（ｒ０）＝φ２（ｒ０）が成立するようにするためには、調整定数Ｃを下記のように設定すればよい。
【数１８】
2018049583.tif 000020
【００５４】
図５は、相互作用ポテンシャルを表す関数φ（ｒ）の形状を示すグラフである。式（１５）のシフティドレナードジョーンズポテンシャルを表す関数φ（ｒ）の形状を曲線ａ０で示し、式（１４）のｂ＝２とした変換則により変換した後のポテンシャルを表す関数φ（ｒ）の形状を曲線ａ１で示す。参考のために、式（１２）のｂ＝２とした変換則により変換した後のポテンシャルを表す関数φ（ｒ）の形状を破線ａ２で示す。
【００５５】
引力領域においては、式（１４）の変換則を適用することにより、ポテンシャルの深さが浅くなり、ポテンシャルの幅が狭くなっている。このように、ポテンシャルの深さを規定するフィッティングパラメータ（式（１１）のε２）を小さくし、ポテンシャルの幅を規定するフィッティングパラメータ（式（１１）のσ２）を小さくして、式（１１）の右辺第２項を不変にすることにより、音速及びヤング率が不変になるような変換が行われる。
【００５６】
斥力領域においては、式（１４）の変換則を適用することにより、変換前のポテンシャルの形状を表す曲線ａ０をポテンシャルの値が小さくなる方向に平行移動することにより、変換後のポテンシャルの形状を表す曲線ａ１が得られる。
【００５７】
図６は、シミュレーションにより求めた動粘度と定数ｂとの関係を示すグラフである。横軸は定数ｂを表し、縦軸は動粘度を単位「ｍ２／ｓ」で表す。定数ｂ＝１のとき、すなわちフィッティングパラメータε及びσの変換前において、シミュレーションの条件を下記の通りとした。
フィッティングパラメータε＝４０４．５Ｋ
フィッティングパラメータσ＝２．６４１×１０−１０ｍ
平衡時の粒子間の距離ｒ０＝２．９６４×１０−１０ｍ
粒子の質量ｍ＝１８．０２Ｄａ
粒子系の温度Ｔ＝３００．０Ｋ
粒子系の圧力Ｐ＝２０．０ａｔｍ（約２．０ＭＰａ）
密度ρは、式（１４）の定数ｂを変更するたびに、Ｔ＝３００．０Ｋ、Ｐ＝２０．０ａｔｍになるように調整した。
【００５８】
図６に示したシミュレーション結果からわかるように、定数ｂを大きくすると動粘度が低下する。例えば、定数ｂを１２８程度にすると、動粘度がフィッティングパラメータ変換前の値に比べて１／３０程度まで低下している。
【００５９】
次に、分子動力学法によるシミュレーションで適用する時間刻み幅について説明する。
一般的に、レナードジョーンズポテンシャルを用いてシミュレーションを行うときの時間刻み幅ｄｔは、粒子の質量をｍとして、以下の式を満たすように設定される。
【数１９】
2018049583.tif 000021
【００６０】
シフティドレナードジョーンズポテンシャルにおいても同様に考えると、時間刻み幅ｄｔは、以下の式を満たすように設定することが好ましい。
【数２０】
2018049583.tif 000022
【００６１】
式（１２）の変換則を用いて相互作用ポテンシャルを表す関数φ（ｒ）を変換した場合、斥力領域におけるフィッティングパラメータの変換を考慮すると、式（２０）は、以下のように変換される。
【数２１】
2018049583.tif 000023
【００６２】
動粘度を低下させるために定数ｂを大きくすると、式（２０）で設定される時間刻み幅ｄｔが小さくなる。これは、式（１２）の変換則を用いてシミュレーションを行う場合に、変換前のフィッティングパラメータを用いたシミュレーションと同等の精度を得るためには、時間刻み幅ｄｔを小さくしなければならないことを意味する。
【００６３】
これに対し、式（１４）の変換則では、斥力領域においてフィッティングパラメータは変換の前後で同一である。このため、式（１４）の変換則を用いてシミュレーションを行う場合には、変換前のフィッティングパラメータを用いたシミュレーションと同等の精度を得るために時間刻み幅ｄｔを小さくする必要はない。
【００６４】
このように、式（１４）の変換則を用いてシミュレーションを行う場合には、定数ｂを大きくして動粘度を低下させるとともに、シミュレーションを行う時間刻み幅ｄｔとして、フィッティングパラメータの変換前のシミュレーションで適用する時間刻み幅をそのまま適用することができる。
【００６５】
次に、図７を参照して実施例によるシミュレーション装置について説明する。
図７は、実施例によるシミュレーション装置の機能ブロック図である。シミュレーション装置としてシミュレーションプログラムを実行するコンピュータを用いることができる。このシミュレーション装置は、中央処理装置（ＣＰＵ）１０、メインメモリ１１、補助記憶装置１２、入力装置１３、及び出力装置１４を含む。補助記憶装置１２は、シミュレーション用のコンピュータプログラムが格納されている記憶媒体を含む。この記憶媒体は、補助記憶装置１２に内蔵されたものでもよく、補助記憶装置１２に着脱可能なリムーバブル媒体でもよい。このコンピュータプログラムがメインメモリ１１にロードされ、ＣＰＵ１０によって実行される。
【００６６】
入力装置１３からシミュレーションに必要な情報が入力される。例えば、ステップＳ１（図１）の相互作用ポテンシャルを表す関数φ（ｒ）を定義する情報、ステップＳ２（図１）で適用されるフィッティングパラメータε及びσの変換則を規定する情報、ステップＳ３（図１）の分子動力学法によるシミュレーションの初期条件、境界条件、シミュレーション対象の粒子系を定義する情報等が入力される。ＣＰＵ１０はシミュレーション結果を出力装置１４に出力する。
【００６７】
上述の各実施例は例示であり、異なる実施例で示した構成の部分的な置換または組み合わせが可能であることは言うまでもない。複数の実施例の同様の構成による同様の作用効果については実施例ごとには逐次言及しない。さらに、本発明は上述の実施例に制限されるものではない。例えば、種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。
【符号の説明】
【００６８】
１０  中央処理装置（ＣＰＵ）
１１  メインメモリ
１２  外部記憶装置
１３  入力装置
１４  出力装置
【図１】
2018049583.tif 000026
【図２】
2018049583.tif 000027
【図３】
2018049583.tif 000028
【図４】
2018049583.tif 000029
【図５】
2018049583.tif 000030
【図６】
2018049583.tif 000031
【図７】
2018049583.tif 000032