(57)【要約】【課題】異方向性を持つ粒子の液相と気相とが混在する状態における粒子の挙動を精度よくシミュレーションする方法を提供する。【解決手段】粒子間の相互作用ポテンシャルが異方向性を持つ粒子の挙動を解析するにあたり、異方向性を持たない粒子間の相互作用ポテンシャルに、異方向性を持たない修正項を付加することにより、フィッティングパラメータを含む修正相互作用ポテンシャルを定義する。修正相互作用ポテンシャルを用いて液相と気相とを含む粒子系の挙動を解析する。解析結果に基づいて得られる飽和蒸気線と、シミュレーション対象の物質の物性値としての飽和蒸気線とを比較して、フィッティングパラメータの値を決定する。決定されたフィッティングパラメータの値を含む修正相互作用ポテンシャルを用いて粒子の挙動を解析する。【選択図】図１
2019164610.tif 000002
【特許請求の範囲】
【請求項１】
粒子間の相互作用ポテンシャルが異方向性を持つ粒子の挙動を解析するシミュレーション方法であって、
異方向性を持たない粒子間の相互作用ポテンシャルに、異方向性を持たない修正項を付加することにより、フィッティングパラメータを含む修正相互作用ポテンシャルを定義し、
前記修正相互作用ポテンシャルを用いて液相と気相とを含む粒子系の挙動を解析し、
解析結果に基づいて得られる飽和蒸気線と、シミュレーション対象の物質の物性値としての飽和蒸気線とを比較して、前記フィッティングパラメータの値を決定し、
決定された前記フィッティングパラメータの値を含む前記修正相互作用ポテンシャルを用いて粒子の挙動を解析するシミュレーション方法。
【請求項２】
前記修正項は、着目する粒子からの距離を変数とするガウシアン関数である請求項１に記載のシミュレーション方法。
【請求項３】
前記相互作用ポテンシャルは、レナードジョーンズ型のポテンシャルであり、前記フィッティングパラメータは、レナードジョーンズ型ポテンシャルのエネルギの次元を持つパラメータと、距離の次元を持つパラメータとを含む請求項１または２に記載のシミュレーション方法。
【請求項４】
さらに、解析結果に基づいて得られる飽和蒸気線と、シミュレーション対象の物質の物性値としての飽和蒸気線とを比較して、シミュレーション対象の粒子の運動方程式に現れる質量の値を決定する請求項１乃至３のいずれか１項に記載のシミュレーション方法。
【請求項５】
圧力条件を含むシミュレーション条件を取得するシミュレーション条件取得部と、
取得された前記圧力条件に応じた粒子間の修正相互作用ポテンシャル、及び運動方程式を用いて粒子の挙動を解析する演算部と、
前記演算部による解析結果を出力する出力制御部と
を有するシミュレーション装置。
【請求項６】
前記圧力条件に応じた粒子間の修正相互作用ポテンシャル及び前記圧力条件に応じた運動方程式を用いた解析結果から求まるシミュレーション対象の系の飽和蒸気線を示す値と、前記圧力条件におけるシミュレーション対象の物質の飽和蒸気線の物性値との誤差が前記物性値の１５％以下である請求項５に記載のシミュレーション装置。
【請求項７】
圧力条件を含むシミュレーション条件を取得する機能と、
取得された前記圧力条件に応じた粒子間の修正相互作用ポテンシャル、及び運動方程式を用いて粒子の挙動を解析する機能と、
解析結果を出力する機能と
をコンピュータに実行させるプログラム。
【請求項８】
前記圧力条件に応じた粒子間の修正相互作用ポテンシャル及び前記圧力条件に応じた運動方程式を用いた解析結果から求まるシミュレーション対象の系の飽和蒸気線を示す値と、前記圧力条件におけるシミュレーション対象の物質の飽和蒸気線の物性値との誤差が前記物性値の１５％以下である請求項７に記載のプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、シミュレーション方法、シミュレーション装置、及びプログラムに関する。
【背景技術】
【０００２】
分子間（粒子間）の相互作用ポテンシャルとしてレナードジョーンズ型ポテンシャルを用いた分子動力学法（ＭＤ法）によるシミュレーションは、不活性ガスの挙動を精度よく再現する。ところが、粒子間の相互作用が異方向性を持つ水分子等の挙動を精度よく再現できない場合がある。水分子等の異方向性を考慮した相互作用ポテンシャルを用いてシミュレーションするには分子の回転を計算する必要があり、計算負荷が著しく増大してしまう。
【０００３】
レナードジョーンズ型ポテンシャルに、静電相互作用ポテンシャルを組み合わせたポテンシャルを作成し、高分子の分子運動シミュレーションに適用可能で、新規材料の評価に好適な分子間相互作用計算手法が提案されている（特許文献１参照）。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特開２０１１−２４８５４１号公報
【発明の概要】
【発明が解決しようとする課題】
【０００５】
ＭＤ法において、従来の相互作用ポテンシャルは、異方向性を持つ分子の相転移（液相から気相、及び気相から液相）を伴う系、例えば蒸気タービン低圧段領域等の水蒸気の挙動を精度よくシミュレーションする目的には適さない。本発明の目的は、異方向性を持つ粒子の液相と気相とが混在する状態における粒子の挙動を精度よくシミュレーションする方法、装置、及びプログラムを提供することである。
【課題を解決するための手段】
【０００６】
本発明の一観点によると、
粒子間の相互作用ポテンシャルが異方向性を持つ粒子の挙動を解析するシミュレーション方法であって、
異方向性を持たない粒子間の相互作用ポテンシャルに、異方向性を持たない修正項を付加することにより、フィッティングパラメータを含む修正相互作用ポテンシャルを定義し、
前記修正相互作用ポテンシャルを用いて液相と気相とを含む粒子系の挙動を解析し、
解析結果に基づいて得られる飽和蒸気線と、シミュレーション対象の物質の物性値としての飽和蒸気線とを比較して、前記フィッティングパラメータの値を決定し、
決定された前記フィッティングパラメータの値を含む前記修正相互作用ポテンシャルを用いて粒子の挙動を解析するシミュレーション方法が提供される。
【０００７】
本発明の他の観点によると、
圧力条件を含むシミュレーション条件を取得するシミュレーション条件取得部と、
取得された前記圧力条件に応じた粒子間の修正相互作用ポテンシャル、及び運動方程式を用いて粒子の挙動を解析する演算部と、
前記演算部による解析結果を出力する出力制御部と
を有するシミュレーション装置が提供される。
【０００８】
本発明のさらに他の観点によると、
圧力条件を含むシミュレーション条件を取得する機能と、
取得された前記圧力条件に応じた粒子間の修正相互作用ポテンシャル、及び運動方程式を用いて粒子の挙動を解析する機能と、
解析結果を出力する機能と
をコンピュータに実行させるプログラムが提供される。
【発明の効果】
【０００９】
異方向性を持つ粒子の液相と気相とが混在する状態における粒子の挙動を精度よくシミュレーションすることができる。
【図面の簡単な説明】
【００１０】
【図１】図１は、修正相互作用ポテンシャルを導出する方法を示すフローチャートである。
【図２】図２は、相互作用ポテンシャルＵ０（ｒ）及び修正相互作用ポテンシャルＵ（ｒ）の一例を示すグラフである。
【図３】図３Ａは、粒子の挙動をＭＤ法により解析した結果の粒子の分布を示す模式図であり、図３Ｂは、相互作用ポテンシャルとしてレナードジョーンズ型ポテンシャルを用いて解析した場合の密度分布の一例を示すグラフである。
【図４】図４Ａは、無次元化した修正相互作用ポテンシャルＵ’を用いて解析を行った結果から求めた無次元飽和蒸気線をプロットしたグラフであり、図４Ｂは、ε＝４．８０×１０−２１Ｊ、σ＝３．１０×１０−１０ｍ、ｍ＝３．１０×１０−２６ｋｇに設定した時の飽和蒸気線を、水の物性値としての飽和蒸気線と比較して示すグラフである。
【図５】図５は、圧力と、フィッティングパラメータε、σ、ｍの好ましい値との対応関係の一例を示すグラフである。
【図６】図６は、実施例によるシミュレーション装置のブロック図である。
【図７】図７は、実施例によるシミュレーション方法のフローチャートである。
【発明を実施するための形態】
【００１１】
本発明の実施例では、相互作用ポテンシャルが異方向性を持つ粒子の挙動のシミュレーションを、異方向性を持たない修正相互作用ポテンシャルを用いて行う。図１〜図４Ｂを参照して、修正相互作用ポテンシャルを導出する方法について説明する。
【００１２】
図１は、修正相互作用ポテンシャルを導出する方法を示すフローチャートである。まず、ステップＳＡ１において、フィッティングパラメータを含む修正相互作用ポテンシャルＵ（ｒ）を、異方向性を持たない粒子間に適用される相互作用ポテンシャルＵ０（ｒ）に、修正項Ｕｃ（ｒ）を付加して以下のように定義する。ここで、ｒは粒子間距離を表す。
【数１】
2019164610.tif 000003
【００１３】
異方向性を持たない粒子間に適用される相互作用ポテンシャルＵ０（ｒ）として、通常、以下に示すレナードジョーンズ型ポテンシャルを用いる。
【数２】
2019164610.tif 000004
ここで、εはエネルギの次元を持つフィッティングパラメータであり、σは長さの次元を持つフィッティングパラメータである。
【００１４】
修正項Ｕｃ（ｒ）として、例えば着目する粒子からの距離ｒを変数とする下記のガウシアン関数を採用する。
【数３】
2019164610.tif 000005
ここで、α、β、γはフィッティングパラメータである。これらのフィッティングパラメータの値は、例えば、異方向性を持つ粒子の相互作用ポテンシャルを全ての方向について平均化したポテンシャル形状を参考にして決定するとよい。例えば、分子動力学法を用いて水分子の挙動をシミュレーションする場合、α＝−１．０、β＝１．１２、γ＝０．２に設定するとよい。
【００１５】
図２は、相互作用ポテンシャルＵ０（ｒ）及び修正相互作用ポテンシャルＵ（ｒ）の一例を示すグラフである。図２の横軸は、粒子間の距離ｒを無次元化した無次元量ｒ／σを表し、縦軸は、ポテンシャルの大きさを無次元化した無次元量Ｕ／εを表す。図２の破線は、無次元化したレナードジョーンズ型ポテンシャルＵｏを示し、実線は、無次元化した修正相互作用ポテンシャルＵを示す。
【００１６】
フィッティングパラメータεは、レナードジョーンズ型の相互作用ポテンシャルＵ０（ｒ）の深さを表し、フィッティングパラメータσは、相互作用ポテンシャルＵ０（ｒ）の大きさが０になる粒子間距離を表す。フィッティングパラメータαは、修正相互作用ポテンシャルＵ（ｒ）の深さに係わるものである。フィッティングパラメータβは、修正相互作用ポテンシャルＵ（ｒ）が極小値をとる距離ｒの値に係わるものである。フィッティングパラメータγは、修正相互作用ポテンシャルＵ（ｒ）の広がりに係わるものである。この段階で、フィッティングパラメータε及びσが未定である。
【００１７】
修正相互作用ポテンシャルＵ（ｒ）を定義した後、ステップＳＡ２（図１）において、修正相互作用ポテンシャルＵ（ｒ）を用いて液相及び気相を含む粒子系の挙動を解析する。以下、この解析方法について説明する。
【００１８】
まず、解析対象の仮想領域内に、粒子の初期条件として粒子密度の高い領域と低い領域とを設定する。解析対象の仮想領域全体について温度制御を行いながらＭＤ法を用いて粒子の挙動を解析する。温度制御には、例えばＮｏｓｅ−Ｈｏｏｖｅｒ法を用いることができる。粒子系の温度が、液相と気相とが混在する温度である場合、仮想領域内に、液相の領域と気相の領域とが現れる。
【００１９】
図３Ａは、粒子の挙動をＭＤ法により解析した結果の粒子の分布を示す模式図である。高さ方向の寸法が最も大きな直方体の仮想領域１０内の高さ方向のほぼ中央部に、粒子密度の高い液相領域１１が現れ、その上下に粒子密度の低い気相領域１２が現れている。
【００２０】
気液界面に対して垂直方向の密度分布は、気液界面に対する垂直方向の座標をｙとして以下の式にフィットさせることができる。
【数４】
2019164610.tif 000006
ここで、ρｇ、ρｌ、ｙｃ、λは、それぞれ気相領域１２の密度、液相領域１１の密度、境界の位置、及び境界の厚さである。
【００２１】
図３Ｂは、相互作用ポテンシャルとしてレナードジョーンズ型ポテンシャルを用いて解析した場合の密度分布の一例を示すグラフである。横軸は、高さ方向の位置を無次元量で表し、縦軸は密度を無次元量で表す。グラフ中の黒丸記号は、解析結果から求まる密度を示し、曲線は近似線を示す。この解析を、いくつかの温度で繰り返し実行することにより、飽和蒸気線の圧力、温度、及び密度を得ることができる。
【００２２】
次に、解析の具体例について説明する。式（１）に示した修正相互作用ポテンシャルＵ（ｒ）を用いてＭＤ法を適用することにより、粒子の挙動を解析する。ＭＤ法においては、以下の運動方程式を解く。
【数５】
2019164610.tif 000007
ここで、ｍは粒子質量、ベクトルｒｉは粒子ｉの位置、ｒｉｊは、粒子ｉと粒子ｊとの距離、ｒｉｊハットは粒子ｉから粒子ｊに向かう単位ベクトルである。
【００２３】
修正相互作用ポテンシャルＵ（ｒ）を用いて式（５）の運動方程式を解くためには、５個の修正相互作用ポテンシャルＵ（ｒ）の形状と定義するフィッティングパラメータα、β、γ、ε、σ、及び運動方程式を定義するフィッティングパラメータｍの値を決定する必要がある。この６個のフィッティングパラメータを同時に調整することは困難であるため、飽和蒸気線の計算を無次元量で行うことで、フィッティングパラメータε、σ、及びｍを計算から除外する。フィッティングパラメータα、β、γの値は、異方向性を持つ相互作用ポテンシャルに基づいて決定しておくとよい。式（１）〜（３）、（５）のエネルギをεで無次元化し、距離をσで無次元化し、質量をｍで無次元化する。無次元化したエネルギをＵ’、無次元化した距離をｒ’で表すと、これらは以下の式で定義される。
【数６】
2019164610.tif 000008
【００２４】
運動方程式及び修正相互作用ポテンシャルは、以下のように書き換えられる。
【数７】
2019164610.tif 000009
【００２５】
このとき、温度、圧力、密度の無次元量Ｔ’、ｐ’、ρ’は、ボルツマン定数をｋＢとして、それぞれ以下の式で表される。
【数８】
2019164610.tif 000010
【００２６】
図４Ａは、無次元化した修正相互作用ポテンシャルＵ’（ｒ’）を用いて解析を行った結果から求めた無次元飽和蒸気線をプロットしたグラフである。図４Ａの横軸は密度の無次元量ρ’を表し、縦軸は圧力の無次元量ｐ’を表す。図４Ａでは、温度の無次元量Ｔ’が１．２０、１．３０、１．４０及び１．５０のときの解析結果から求まった無次元飽和蒸気線上の点を示している。
【００２７】
フィッティングパラメータε、σ、ｍに種々の値を設定してスケーリングすることにより、次元を持つ飽和蒸気線を描くことができる。
【００２８】
図４Ｂは、ε＝４．８０×１０−２１Ｊ、σ＝３．１０×１０−１０ｍ、ｍ＝３．１０×１０−２６ｋｇとしたときの解析結果から求まる飽和蒸気線を、水の物性値としての飽和蒸気線と比較して示すグラフである。横軸は密度ρを単位「ｋｇ／ｍ３」で表し、縦軸は圧力ｐを単位「ｋＰａ」で表す。グラフ中の黒四角記号は、特定の温度のときの解析結果から求めた飽和蒸気線上の値を示し、黒丸記号は、水の物性値を示す。
【００２９】
図４Ｂに示した例では、圧力が１０００ｋＰａ以下の領域において、計算結果から求まる飽和蒸気線と物性値としての飽和蒸気線とがよく一致している。例えば、解析結果から得られた飽和蒸気線を示す値と、物性値との誤差が、物性値の１５％以下である。ただし、圧力が１０００ｋＰａより高い領域では、解析結果から求めた飽和蒸気線の値と物性値との誤差が大きい。
【００３０】
フィッティングパラメータε、σ、ｍに上記の値を設定した粒子系を用いて、圧力が１０００ｋＰａ以下の範囲において、液相と気相とが混在する水分子の挙動を精度よくシミュレーションすることができることがわかる。ところが、圧力がそれよりも高い範囲では、シミュレーションの精度が低下する。
【００３１】
フィッティングパラメータε、σ、ｍの値を変化させると、圧力が１０００ｋＰａよりも高い範囲で、解析結果から求まる飽和蒸気線と、物性値としての飽和蒸気線との誤差を小さくすることが可能である。シミュレーションする系の圧力条件に応じて、解析結果から求まる飽和蒸気線と、物性値としての飽和蒸気線との誤差が小さくなるように、フィッティングパラメータε、σ、ｍの値を設定するとよい。
【００３２】
圧力に応じたフィッティングパラメータε、σ、ｍの好ましい値は、予め種々のシミュレーションを行って決めておくとよい。
【００３３】
図５は、圧力と、フィッティングパラメータε、σ、ｍの好ましい値との対応関係の一例を示すグラフである。シミュレーションすべき系の圧力の範囲（圧力条件）が決まると、図５から、その圧力条件に適したフィッティングパラメータε、σ、ｍの値を決定し、その値を用いてシミュレーションを行えばよい。
【００３４】
次に、図６を参照して、実施例によるシミュレーション装置について説明する。
図６は、実施例によるシミュレーション装置のブロック図である。本実施例によるシミュレーション装置は、処理装置３０、入力装置３８、及び出力装置３９を含む。処理装置３０は、シミュレーション条件取得部３１、演算部３２、出力制御部３３、及び圧力−パラメータ対応関係記憶部３４を含む。圧力−パラメータ対応関係記憶部３４に、例えば図５に示した圧力とパラメータとの対応関係が記憶されている。
【００３５】
図６に示す各ブロックは、ハードウェア的には、コンピュータの中央処理ユニット（ＣＰＵ）をはじめとする素子や機械装置で実現することができ、ソフトウェア的にはコンピュータプログラム等によって実現することができる。図６では、ハードウェア及びソフトウェアの連携によって実現される機能ブロックが示されている。従って、これらの機能ブロックは、ハードウェア及びソフトウェアの組み合わせによって、種々の態様で実現することが可能である。
【００３６】
処理装置３０は入力装置３８及び出力装置３９と接続される。入力装置３８は、処理装置３０で実行される処理に関係するユーザからのコマンド及びデータの入力を受ける。入力装置３８として、例えばユーザが操作を行うことにより入力を行うキーボードやマウス、インターネット等のネットワークを介して入力を行う通信装置、ＣＤ、ＤＶＤ等の記録媒体から入力を行う読取装置等を用いることができる。
【００３７】
シミュレーション条件取得部３１は、入力装置３８を介してシミュレーション条件を取得する。シミュレーション条件には、シミュレーションに必要な種々の情報が含まれる。例えば、シミュレーション対象の粒子（分子）に関する物理量の初期条件、境界条件、シミュレーション対象の圧力条件等が含まれる。
【００３８】
演算部３２は、シミュレーション対象の圧力条件、及び圧力−パラメータ対応関係記憶部３４に記憶されている対応関係に基づいて、ＭＤ法の実行時に用いる修正相互作用ポテンシャルＵ（ｒ）のフィッティングパラメータε、σ、及び運動方程式のフィッティングパラメータｍの値を決定する。さらに、演算部３２は、決定されたフィッティングパラメータε、σ、ｍの値を用いてＭＤ法により粒子の挙動をシミュレーションする。
【００３９】
出力制御部３３は、シミュレーション結果を出力装置３９に出力する。例えば、粒子の位置の変動を、出力装置３９の表示画面に図形で表示する。
【００４０】
次に、図７を参照して、実施例によるシミュレーション方法について説明する。
図７は、実施例によるシミュレーション方法のフローチャートである。ステップＳＢ１において、シミュレーション条件取得部３１（図６）がシミュレーション条件、及びシミュレーション対象の系の圧力条件を取得する。
【００４１】
次に、ステップＳＢ２において、演算部３２が、取得された圧力条件、及び圧力−パラメータ対応関係記憶部３４（図６）に記憶されている対応関係から、フィッティングパラメータε、σ、ｍの値を決定する。さらに、ステップＳＢ３において、演算部３２が、決定されたフィッティングパラメータε、σ、ｍの値を用いてＭＤ法により粒子の挙動をシミュレーションする。
【００４２】
その後、ステップＳＢ４において、出力制御部３３がシミュレーション結果を出力装置３９に出力する。
【００４３】
次に、上記実施例の持つ優れた効果について説明する。
実施例による方法では、シミュレーション結果から求まる飽和蒸気線が、シミュレーション対象の物性値としての飽和蒸気線とほぼ一致する修正相互作用ポテンシャルＵ（ｒ）をＭＤ法に適用している。このため、液相と気相とが混在する系のシミュレーション精度を高めることができる。従来のレナードジョーンズポテンシャルをＭＤ法に適用した場合では、フィッティングパラメータをどのように設定しても、シミュレーション結果から求まる飽和蒸気線と、物性値としての飽和蒸気線との誤差を小さくすることができなかった。実施例による修正相互作用ポテンシャルＵ（ｒ）を用いることにより、シミュレーション結果から求まる飽和蒸気線と、物性値としての飽和蒸気線との誤差を、物性値の１５％以下にすることが可能になる。
【００４４】
実施例では、修正相互作用ポテンシャルＵ（ｒ）として、異方向性を持たないものを用いるため、異方向性を持つ相互作用ポテンシャルを用いてシミュレーションを行う場合に比べて、計算負荷を軽減させることができる。
【００４５】
シミュレーションに用いるフィッティングパラメータε、σ、ｍの値を変えることにより、種々の圧力に対応したシミュレーションを行うことが可能である。例えば、圧力が１０００ｋＰａ以下の領域に対応するフィッティングパラメータε、σ、ｍを用いることにより、蒸気タービンの低圧段領域での水蒸気及び水の挙動を精度よくシミュレーションすることができる。
【００４６】
次に、上記実施例の変形例について説明する。
上記実施例では、シミュレーション対象として水分子を取り上げたが、上記実施例は、その他の異方向性を持つ分子の挙動のシミュレーションにも適用することが可能である。
【００４７】
上記実施例では、式（３）のフィッティングパラメータα、β、γの値を固定したが、フィッティングパラメータα、β、γについても種々の値に設定して予め解析を行い、物性値としての飽和蒸気線をよく再現する条件を見出しておくとよい。
【００４８】
上記実施例では、修正相互作用ポテンシャルＵ（ｒ）をＭＤ法に適用したが、修正相互作用ポテンシャルＵ（ｒ）を、くりこみ群分子動力学法（ＲＭＤ法）に適用してもよい。この場合には、ＭＤ法を用いるときのフィッティングパラメータε、σ、ｍの値を、これらのパラメータのくりこみ変換則に従って変換すればよい。
【００４９】
上述の各実施例は例示であり、異なる実施例で示した構成の部分的な置換または組み合わせが可能であることは言うまでもない。複数の実施例の同様の構成による同様の作用効果については実施例ごとには逐次言及しない。さらに、本発明は上述の実施例に制限されるものではない。例えば、種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。
【符号の説明】
【００５０】
１０  シミュレーション対象領域
１１  液相領域
１２  気相領域
３０  処理装置
３１  シミュレーション条件取得部
３２  演算部
３３  出力制御部
３４  圧力−パラメータ対応関係記憶部
３８  入力装置
３９  出力装置
【図１】
2019164610.tif 000011
【図２】
2019164610.tif 000012
【図３】
2019164610.tif 000013
【図４】
2019164610.tif 000014
【図５】
2019164610.tif 000015
【図６】
2019164610.tif 000016
【図７】
2019164610.tif 000017