JP7086483B2 - シミュレーション方法、シミュレーション装置及びプログラム - Google Patents

Description

本発明は、分子動力学法またはくりこみ群分子動力学法を用いたシミュレーション方法、シミュレーション装置及びプログラムに関する。

蒸気タービンの低圧段で蒸気の温度が低下すると凝縮が起こり、水滴が形成される。この水滴が動翼に衝突して、動翼にエロージョンが生じる。このエロージョンを抑制するために、蒸気タービン内の蒸気及び水滴の挙動を知ることが重要である。従来、蒸気等の流体の流れ場のシミュレーション解析では、流体を連続体として扱っていた（例えば、特許文献１）。このようなシミュレーション解析では、気体から液体への相変化を伴う流れ場の詳細な挙動を把握することは困難である。

分子動力学法またはくりこみ群分子動力学によるシミュレーション解析を行うことにより、流体の挙動を解析する手法が提案されている（特許文献２）。

特開２０１７－００４１０３号公報 特開２０１６－２１８７６７号公報

分子動力学法を用いた従来のシミュレーション解析では、流体の流入及び流出を伴うような系を取り扱うことが困難であった。本発明の目的は、流体の流入及び流出を伴う系を、分子動力学法を用いて解析するシミュレーション方法、シミュレーション装置及びプログラムを提供することである。本明細書において、くりこみ群分子動力学法は、広義の意味で分子動力学法ということもできる。本明細書において、分子動力学法及びくりこみ群分子動力学法を単に「分子動力学法」という。

本発明の一観点によると、
流出入界面を持つ流れ場を解析領域とし、流れ場内の流体を複数の粒子の集合体として分子動力学法を用いてシミュレーションする方法であって、
前記流出入界面に、境界領域を介して熱浴を接続し、前記熱浴と前記境界領域との間、及び前記境界領域と前記解析領域との間で粒子の移動を許容した解析モデルにおいて、
前記流出入界面における圧力目標値を取得し、
粒子の状態が時間発展したときに前記境界領域内の圧力が前記圧力目標値に維持されるように、前記熱浴内の圧力を制御するシミュレーション方法が提供される。

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

熱浴内の圧力を制御することにより、解析領域内の流出入界面における圧力を圧力目標値に維持することができる。さらに、解析領域と熱浴との間に圧力制御の対象とならない境界領域を配置することにより、解析領域の流出入界面における圧力を圧力目標値に維持する精度を高めることができる。

図１は、実施例によるシミュレーション装置のブロック図である。 図２は、実施例によるシミュレーション方法で解析する対象となる解析モデルの一例を示す概略図である。 図３は、実施例によるシミュレーション装置の処理部が実行する処理手順のフローチャートである。 図４は、流出入界面における圧力制御の手順を示すフローチャートである。 図５Ａは、第１熱浴内の圧力制御（図４のステップＳ４５）の手順を示すフローチャートであり、図５Ｂは、個数ｄＮが正の場合の粒子の追加の様子を示す第１熱浴の模式図であり、図５Ｃは、個数ｄＮが負の場合の粒子の除去の様子を示す第２熱浴の模式図である。 図６Ａは、実施例及び比較例による方法でシミュレーションを行って求めた圧力のノズル軸方向に関する分布を示すグラフであり、図６Ｂは、図６Ａの一部を拡大したグラフである。 図７は、他の実施例によるシミュレーション方法で取り扱う解析モデルの斜視図である。

図１～図６Ｂを参照して、実施例によるシミュレーション装置、及びシミュレーション方法について説明する。

図１は、実施例によるシミュレーション装置のブロック図である。実施例によるシミュレーション装置は、入力部１０、処理部１１、出力部１２、及び記憶部１３を含む。入力部１０から処理部１１にシミュレーション条件等が入力される。さらに、オペレータから入力部１０に各種指令（コマンド）等が入力される。入力部１０は、例えば通信装置、リムーバブルメディア読取装置、キーボード等で構成される。

処理部１１は、入力されたシミュレーション条件及び指令に基づいて分子動力学法またはくりこみ群分子動力学法（以下、単に分子動力学法という。）を用いたシミュレーションを行う。さらに、シミュレーション結果を出力部１２に出力する。シミュレーション結果には、シミュレーション対象物である粒子系の物理量の時間的変化を表す情報が含まれる。処理部１１は、例えばコンピュータを含み、分子動力学法によるシミュレーションをコンピュータに実行させるためのプログラムが記憶部１３に記憶されている。出力部１２は、通信装置、リムーバブルメディア書込み装置、ディスプレイ等を含む。

図２は、実施例によるシミュレーション方法で解析する対象となる解析モデルの一例を示す概略図である。例えば壁面２３、及び一対の流出入界面２１、２２を持つ四角柱状の解析領域２０が定義される。一方の流出入界面２１から解析領域２０内に流体、例えば水蒸気が流入し、他方の流出入界面２２から外部に流出する流れ場が解析領域２０内に形成される。この流体を複数の粒子の集合体で表し、分子動力学法を用いて粒子の挙動を解析することにより、解析領域２０内の流れ場の解析を行う。境界条件として、一方の流出入界面２１における温度、圧力の目標値（温度目標値Ｔ１ｏ、圧力目標値Ｐ１ｏ）、及び他方の流出入界面２２における温度、圧力の目標値（温度目標値Ｔ２ｏ、圧力目標値Ｐ２ｏ）が与えられる。

流出入界面２１に、境界領域３０を介して第１熱浴３１が接続されている。第１熱浴３１と境界領域３０との間、及び境界領域３０と解析領域２０との間で、粒子５０の移動が許容されている。もう一方の流出入界面２２に第２熱浴４０が接続されている。第２熱浴４０と解析領域２０との間で粒子５０の移動が許容されている。

流出入界面２１、２２に接する面を通って解析領域２０に粒子５０が流入し、または解析領域２０から粒子５０が流出する。第１熱浴３１の表面のうち境界領域３０との界面以外の表面、境界領域３０の表面のうち、第１熱浴３１との境界及び流出入界面２１以外の表面、第２熱浴４０の表面のうち流出入界面２２以外の表面には、反射境界条件を適用する。反射境界に接触した粒子５０は、接触時とは反対の面法線方向の速度成分をもって反射される。解析領域２０の壁面２３には、例えば周期境界条件または反射境界条件を適用する。第１熱浴３１及び境界領域３０内の圧力を第２熱浴４０内の圧力より高く設定すると、解析領域２０内に流出入界面２１から２２に向かう流れ場が形成される。

図２に示した解析モデルを通常の分子動力学法により解析すると、流れ場の上流側の流出入界面２１に接続された境界領域３０及び第１熱浴３１内の粒子５０は、解析領域２０内に流出することにより、時間の経過とともに減少する。逆に、流れ場の下流側の流出入界面２２に接続された第２熱浴４０内の粒子５０は、解析領域２０から流入することにより、時間の経過とともに増加する。粒子５０の増減により、境界領域３０及び第１熱浴３１内の圧力は時間の経過とともに低下し、第２熱浴４０内の圧力は時間の経過とともに上昇する。このため、流出入界面２１及び流出入界面２２における圧力を一定に維持することはできない。以下に説明する実施例では、流出入界面２１及び流出入界面２２における圧力を一定の圧力目標値Ｐ１ｏ、Ｐ２ｏに維持する処理を実行する。

図３は、実施例によるシミュレーション装置の処理部１１（図１）が実行する処理手順のフローチャートである。

まず、処理部１１は、入力部１０に入力されたシミュレーションの初期条件、境界条件等のシミュレーション条件を取得する（ステップＳ１）。境界条件には、解析領域２０、境界領域３０、第１熱浴３１、及び第２熱浴４０の形状及び大きさ、流出入界面２１、２２における圧力目標値Ｐ１ｏ、Ｐ２ｏ、及び温度目標値Ｔ１ｏ、Ｔ２ｏが含まれる。初期条件には、粒子５０の位置及び速度を示す情報が含まれる。さらに、シミュレーション条件として、粒子５０の質量と大きさ、粒子間の相互作用ポテンシャルを規定する情報、時間発展させる時間刻み幅等の情報が含まれる。粒子間の相互作用ポテンシャルとして、例えばレナードジョーンズポテンシャルを適用することができる。相互作用ポテンシャルを規定する情報には、例えばレナードジョーンズポテンシャルのフィッティングパラメータが含まれる。

シミュレーション条件を取得すると、処理部１１は初期条件に基づいて複数の粒子５０を第１熱浴３１、境界領域３０、解析領域２０、及び第２熱浴４０内に配置する。粒子５０の間の相互作用ポテンシャルに基づいて運動方程式を解くことにより、粒子５０の次状態を算出する（ステップＳ３）。具体的には、１タイムステップ後における粒子５０の位置及び速度を算出する。

粒子５０の次状態が求まると、流出入界面２１、２２における圧力Ｐ１、Ｐ２を圧力目標値に維持するための圧力制御を行う（ステップＳ４）。以下、圧力制御について簡単に説明する。なお、圧力制御の詳細については、後に、図４及び図５Ａ～図５Ｃを参照して説明する。

まず、流出入界面２１における圧力制御について説明する。粒子５０の最新の状態（すなわちステップＳ３で求められた粒子５０の次状態）における境界領域３０（図２）内の圧力Ｐ１ｂ、及び第１熱浴３１内の圧力Ｐ１ｈを計算する。境界領域３０内の圧力Ｐ１ｂと圧力目標値Ｐ１ｏとを比較し、この比較結果、及び第１熱浴３１内の最新の状態における圧力Ｐ１ｈに基づいて、第１熱浴３１の圧力目標値Ｐ１ｈｏを決定する。第１熱浴３１内の圧力が圧力目標値Ｐ１ｈｏになるように、第１熱浴３１内に粒子５０を追加、または第１熱浴３１内から粒子５０を除去する。

次に、流出入界面２２における圧力制御について説明する。ステップＳ３で求められた粒子５０の最新の状態における第２熱浴４０（図２）内の圧力Ｐ２ｈを計算する。第２熱浴４０の圧力Ｐ２ｈと圧力目標値Ｐ２ｏとを比較し、この比較結果に基づいて、第２熱浴４０内の圧力が圧力目標値Ｐ２ｏになるように、第２熱浴４０内に粒子５０を追加、または第２熱浴４０内から粒子５０を除去する。

圧力制御を行った後、流出入界面２１、２２における温度Ｔ１、Ｔ２を、温度目標値に維持するための温度制御を行う（ステップＳ５）。具体的には、境界領域３０及び第１熱浴３１内の粒子５０の温度を流出入界面２１における目標温度に維持し、第２熱浴４０内の粒子５０の温度を流出入界面２２における目標温度に維持する制御を行う。温度制御には、例えば速度スケーリング法を用いることができる。

圧力制御及び温度制御を行った後、タイムステップを更新する（ステップＳ６）。具体的には、ステップＳ３で算出された粒子５０の次状態に対して圧力制御及び温度制御を行った後の粒子５０の状態を、現在の状態として設定する。

解析を終了するまで、ステップＳ３からステップＳ６までの処理を繰り返す（ステップＳ７）。解析を終了する場合には、シミュレーショ結果を出力部１２（図１）に出力する（ステップＳ８）。出力部１２に出力する情報には、流出入界面２１、２２における圧力の時間変化、解析領域２０内の粒子の状態の時間変化に関する情報を含めるとよい。

次に、図４～図５Ｃを参照して、流出入界面２１（図２）における圧力制御（図３のステップＳ４）について説明する。
図４は、流出入界面２１（図２）における圧力制御の手順を示すフローチャートである。まず、境界領域３０内の圧力を計算する（ステップＳ４１）。境界領域３０内の圧力は、例えばビリアル定理を用いて計算することができる。

次に、現時点のタイムステップにおいて圧力制御を実行するか否かを判定する（ステップＳ４２）。例えば、圧力制御は、数百タイムステップごとに１回行う。圧力制御を実行しない場合は、圧力制御の処理を終了して、図３に示したフローチャートに戻る。圧力制御を実行する場合は、第１熱浴３１の圧力目標値Ｐ１ｈｏを更新するか否かを判定する（ステップＳ４３）。圧力目標値Ｐ１ｈｏの更新は、例えば、圧力制御を１０～１００回行うごとに１回行う。

圧力目標値Ｐ１ｈｏを更新しない場合には、第１熱浴３１の現時点の圧力目標値Ｐ１ｈｏに基づいて第１熱浴３１の圧力を制御する（ステップＳ４５）。圧力目標値Ｐ１ｈｏの初期値として、流出入界面２１における圧力目標値Ｐ１ｏを採用する。圧力目標値Ｐ１ｈｏを更新する場合には、境界領域３０の圧力Ｐ１ｂの平均値に基づいて、第１熱浴３１の圧力目標値Ｐ１ｈｏを更新する（ステップＳ４４）。境界領域３０の圧力Ｐ１ｂの平均値は、ステップＳ４１で求められた過去の複数のタイムステップにおける圧力Ｐ１ｂを平均することにより求められる。第１熱浴３１の圧力目標値Ｐ１ｈｏは、境界領域３０の圧力Ｐ１ｂの平均値、流出入界面２１における圧力目標値Ｐ１ｏ、及び第１熱浴３１の現時点の圧力目標値Ｐ１ｈｏに基づいて決定する。例えば、圧力目標値Ｐ１ｈｏは以下の式に基づいて更新する。

式（１）において左辺のＰ１ｈｏは更新後の圧力目標値であり、右辺のＰ１ｈｏは現時点（更新前）の圧力目標値である。

圧力目標値Ｐ１ｈｏを更新した後、第１熱浴３１の更新後の圧力目標値Ｐ１ｈｏに基づいて第１熱浴３１の圧力を制御する（ステップＳ４５）。例えば、境界領域３０の圧力Ｐ１ｂが圧力目標値Ｐ１ｏより低ければ、その程度に応じて、第１熱浴３１の圧力目標値Ｐ１ｈｏを、最新の状態における第１熱浴３１の圧力Ｐ１ｈより高くする。逆に、境界領域３０の圧力Ｐ１ｂが圧力目標値Ｐ１ｏより高ければ、その程度に応じて、第１熱浴３１の圧力目標値Ｐ１ｈｏを、最新の状態における第１熱浴３１の圧力Ｐ１ｈより低くする。このように、第１熱浴３１の圧力目標値Ｐ１ｈｏを一定間隔で修正する。

次に、流出入界面２２（図２）における圧力制御（図３のステップＳ４）について説明する。流出入界面２２においては、第２熱浴４０内の圧力目標値として流出入界面２２における圧力目標値Ｐ２ｏを用い、第２熱浴４０内の圧力を制御する。

次に、図５Ａ～図５Ｃを参照して、第１熱浴３１内の圧力制御について説明する。
図５Ａは、第１熱浴３１内の圧力制御（図４のステップＳ４５）の手順を示すフローチャートである。ステップＳＡ４で算出された圧力（粒子５０の次状態における圧力）と圧力の目標値との比較結果に基づいて、第１熱浴３１及び第２熱浴４０に対して追加または除去する粒子５０の数を算出する（ステップＳ４５１）。

ステップＳ３で算出した最新の状態における第１熱浴３１内の粒子５０の総数をＮとしたとき、追加すべき粒子の個数ｄＮを、以下の式で算出する。

ここで、ｆｌｏｏｒ関数は、小数点以下を切り捨てて整数にする関数である。なお、小数点以下を切り捨てる代わりに切り上げてもよいし、小数点第１位を四捨五入して整数にしてもよい。

算出された個数ｄＮに相当する数の粒子を、第１熱浴３１に対して追加または除去する（ステップＳ４５２）。具体的には、個数ｄＮが正のとき粒子５０を追加し、個数ｄＮが負のとき粒子５０を除去し、個数ｄＮが０のとき粒子５０の追加も除去も行わない。第２熱浴４０に対しても、第１熱浴３１と同様に粒子５０の追加または除去を行う。

図５Ｂは、個数ｄＮが正の場合の粒子５０の追加の様子を示す第１熱浴３１の模式図である。図５Ｂにおいて、直近の粒子５０の位置を破線で表し、最新の状態（ステップＳ３で求められた次状態）の粒子５０の位置を実線で表している。直近の状態から最新の状態に遷移するときに、２個の粒子５０ｃが第１熱浴３１から流出している。最新の状態では、第１熱浴３１内の粒子数が減少するため、圧力が低下する。このため、Ｐ１ｈｏ＞Ｐ１ｈが成立し、個数ｄＮが正になる。一例として、ｄＮ＝２の場合、図５Ｂの下図に示すように、第１熱浴３１内に２つの粒子５０ａを追加する。

新たに追加された粒子５０ａと、既に存在する粒子５０との距離が短すぎると、レナードジョーンズポテンシャルによる大きな斥力が両者に作用する。その結果、次のタイムステップで粒子が急激に加速され、計算が破綻する可能性が生じる。新たに追加した粒子５０ａと、既に存在する粒子５０との間に作用する斥力が所定の許容上限値を超える場合には、新たに追加した粒子５０ａを再配置する。

図５Ｃは、個数ｄＮが負の場合の粒子５０の除去の様子を示す第２熱浴４０の模式図である。図５Ｃにおいて、直近の状態における粒子５０の位置を破線で表し、最新の状態の粒子５０の位置を実線で表している。直近の状態から最新の状態に遷移するときに、２個の粒子５０ｄが第２熱浴４０に流入している。最新の状態では、第２熱浴４０内の粒子数が増加するため、圧力が上昇する。このため、Ｐ２ｈｏ＜Ｐ２ｈが成立し、個数ｄＮが負になる。一例として、ｄＮ＝－２の場合、図５Ｃの下図に示すように、第２熱浴４０内から２つの粒子５０ｂを除去する。

次に、上記実施例の優れた効果について説明する。
本願発明者による種々のシミュレーション実験によると、圧力を調整するために粒子５０の追加を行っている第１熱浴３１と、解析領域２０とを直接接続すると、両者の界面において圧力が不連続に変化する現象が発生し得ることがわかった。両者の界面で圧力が不連続に変化すると、第１熱浴３１内の圧力が圧力目標値Ｐ１ｏに維持されていても、流出入界面２１から見て解析領域２０側の微小な領域内の圧力が圧力目標値Ｐ１ｏからずれてしまう。

上記実施例では、境界領域３０に対しては粒子５０の追加及び除去を行わない。このため、境界領域３０と解析領域２０との間の流出入界面２１において圧力の不連続な変化は生じない。境界領域３０内の圧力Ｐ１ｂは、流出入界面２１における圧力目標値Ｐ１ｏに維持される。このため、粒子５０の追加を行っている第１熱浴３１と、粒子の追加を行わない境界領域３０との界面において圧力が不連続に変化したとしても、流出入界面２１から見て解析領域２０側の微小な領域内の圧力を圧力目標値Ｐ１ｏに維持することができる。

上述のように、上記実施例は、圧力境界条件を伴う非平衡分子動力学によるシミュレーションに適用することにより、圧力境界条件をシミュレーションに高精度に反映させることができる。

次に、図６Ａ及び図６Ｂを参照して、実施例の優れた効果を確認するために行ったシミュレーション及びその結果について説明する。

実施例による方法及び比較例による方法で、ラバルノズル（収束拡大ノズル）内を流れる流体の圧力分布をシミュレーションにより求めた。実施例による方法では、ラバルノズルの入口に、図２に示した第１熱浴３１と解析領域２０を接続した。比較例による方法では、ラバルノズルの入口に第１熱浴３１を直接接続した。ラバルノズルの出口は開放した。

図６Ａは、実施例及び比較例による方法でシミュレーションを行って求めた圧力のノズル軸方向に関する分布を示すグラフである。図６Ｂは、図６Ａの一部を拡大したグラフである。横軸は、第１熱浴３１の端面からラバルノズルの出口までの長さで正規化した距離を表し、縦軸は入口における圧力目標値Ｐ１ｏで正規化した圧力を表す。横軸の左端が第１熱浴３１の端面に対応し、右端がラバルノズルの出口に対応する。中空の丸記号及び中実の丸記号が、それぞれ実施例及び比較例によるシミュレーション結果を示す。

図６Ａに示したように、入口から出口に向かって圧力が低下している。図６Ｂに示した左端から３個目までの３個の中空の丸記号及び中実の丸記号が、第１熱浴３１内の圧力を示し、左から４個目の中空の丸記号が、境界領域３０内の圧力を示している。左端から５個及びそれよりも右側の中空の丸記号、及び左端から４個目、及びそれよりも右側の中実の丸記号が、解析領域２０内の圧力を示している。

比較例においては、解析領域２０の左端における正規化圧力が１より低くなっている。すなわち、解析領域２０の左端における圧力が圧力目標値Ｐ１ｏに維持されていない。これに対し実施例においては、境界領域３０内の正規化圧力が１になるように第１熱浴３１内の圧力が制御されているため、境界領域３０内の正規化圧力がほぼ１に維持されている。境界領域３０と解析領域２０との境界における圧力の変化は緩やかであるため、解析領域２０の左端における正規化圧力もほぼ１に維持されている。すなわち、解析領域２０の左端における圧力が、圧力目標値Ｐ１ｏに維持されている。

上記シミュレーション実験により、解析領域２０の圧力境界条件が規定されている面の圧力を、圧力目標値に維持することが可能であることが確認された。

次に、上記実施例の変形例について説明する。
上記実施例では、四角柱状の解析領域２０の両端の流出入界面２１、２２（図２）に圧力境界条件を適用しているが、その他の面にも圧力境界条件を適用することが可能である。また、上記実施例では、流れ場の上流側の流出入界面２１に境界領域３０を接続したが、下流側の流出入界面２２に境界領域３０を接続してもよい。

上記実施例では、解析領域２０の形状を四角柱状としているが、その他の任意の形状にすることも可能である。また、流体と相互作用して姿勢を変えたり移動したりする剛体を配置してもよい。

また、上記実施例では、粒子間の相互作用ポテンシャルとして、レナードジョーンズポテンシャルを適用しているが、その他のポテンシャル、例えばモースポテンシャル等を適用してもよい。

圧力目標値Ｐ１ｈｏを更新する処理（図４のステップＳ４４）を実行する周期、及び第１熱浴３１の圧力制御（図４のステップＳ４５）を実行する周期は、任意に設定することができる。これらの処理を実行する周期を短くしすぎると計算負荷が増大する。逆に、これらの処理を実行する周期を長くしすぎると、解析領域２０の圧力境界条件が設定されている流出入界面２１（図２）における圧力と、圧力目標値Ｐ１ｏとのずれが大きくなる。圧力境界条件が設定されている流出入界面２１における圧力と、圧力目標値Ｐ１ｏとのずれが許容範囲内に納まるように、これらの処理を実行する周期を設定することが好ましい。

次に、図７を参照して他の実施例によるシミュレーション方法について説明する。以下、図１～図５に示した実施例によるシミュレーション方法と共通の構成については説明を省略する。図１～図５に示した実施例では、解析領域２０の流出入界面２１における圧力目標値Ｐ１ｏが面内に一定であるが、本実施例では、圧力目標値Ｐ１ｏが面内の場所によって異なっている。

図７は、本実施例によるシミュレーション方法で取り扱う解析モデルの斜視図である。本実施例においても、図１～図５に示した実施例と同様に、解析領域２０に流出入界面２１、２２が画定されている。流出入界面２１が、境界領域３０を介して第１熱浴３１に接続されている。本実施例では、流出入界面２１が複数の区画２１ａに区分されている。流出入界面２１の複数の区画２１ａに対応して、境界領域３０が複数の境界領域セル３０ａに区分され、第１熱浴３１が複数の第１熱浴セル３１ａに区分されている。複数の第１熱浴セル３１ａは、それぞれ対応する境界領域セル３０ａを介して解析領域２０に接続される。

流出入界面２１における圧力目標値Ｐ１ｏの面内分布を反映するように、複数の区画２１ａの各々に圧力目標値Ｐ１ｏが設定される。区画２１ａごとに、当該区画２１ａの圧力目標値Ｐ１ｏが異なっている。隣接する２つの境界領域セル３０ａの界面、隣接する２つの第１熱浴セル３１ａの界面、相互に対応する境界領域セル３０ａと第１熱浴セル３１ａとの界面は、粒子の通過を許容しており、これらの界面を横切って粒子が自由に移動する。

図４のステップＳ４１では、複数の境界領域セル３０ａごとに圧力を計算する。ステップＳ４４では、複数の境界領域セル３０ａのそれぞれについて圧力の平均値を計算する。さらに、第１熱浴セル３１ａのそれぞれについて、対応する境界領域セル３０ａの圧力の平均値、及び対応する区画２１ａにおける圧力目標値Ｐ１ｏに基づいて圧力目標値Ｐ１ｈｏを更新する。ステップＳ４５では、第１熱浴セル３１ａごとに圧力制御を行う。すなわち、図５のステップＳ４５１において、第１熱浴セル３１ａごとに追加または除去する粒子数を算出し、ステップＳ４５２において、第１熱浴セル３１ａごとに粒子の追加または除去を行う。

流出入界面２１を複数の区画２１ａに区分する条件は、シミュレーション条件の一つとして入力部１０（図１）に入力される。入力部１０に入力された条件は、ステップＳ１（図３）で処理部１１（図１）が取得する。処理部１１は、入力されたシミュレーション条件に基づいて、境界領域３０を複数の境界領域セル３０ａに区分し、さらに、第１熱浴３１を複数の第１熱浴セル３１ａに区分する。

次に、第２実施例の優れた効果について説明する。
第２実施例では、圧力境界条件が適用されている流出入界面２１における圧力目標値Ｐ１ｏが面内で一様ではない場合にも、圧力境界条件を満たすようにシミュレーションを行うことができる。

上述の各実施例は例示であり、異なる実施例で示した構成の部分的な置換または組み合わせが可能であることは言うまでもない。複数の実施例の同様の構成による同様の作用効果については実施例ごとには逐次言及しない。さらに、本発明は上述の実施例に制限されるものではない。例えば、種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。

１０ 入力部
１１ 処理部
１２ 出力部
１３ 記憶部
２０ 解析領域
２１ 流出入界面
２１ａ 流出入界面の区画
２２ 流出入界面
２３ 壁面
３０ 境界領域
３０ａ 境界領域セル
３１ 第１熱浴
３１ａ 第１熱浴セル
４０ 第２熱浴
５０、５０ａ、５０ｂ、５０ｃ、５０ｄ 粒子
Ｐ１ 流出入界面における圧力
Ｐ１ｂ 境界領域の圧力
Ｐ１ｈ 第１熱浴の圧力
Ｐ１ｈｏ 第１熱浴の圧力目標値
Ｐ１ｏ 流出入界面における圧力目標値
Ｐ２ 流出入界面における圧力
Ｐ２ｈ 第２熱浴の圧力
Ｐ２ｏ 流出入界面における圧力目標値

Claims (9)

1. 流出入界面を持つ流れ場を解析領域とし、流れ場内の流体を複数の粒子の集合体として分子動力学法を用いてシミュレーションする方法であって、
   前記流出入界面に、境界領域を介して熱浴を接続し、前記熱浴と前記境界領域との間、及び前記境界領域と前記解析領域との間で粒子の移動を許容した解析モデルにおいて、
   前記流出入界面における圧力目標値を取得し、
   粒子の状態が時間発展したときに前記境界領域内の圧力が前記圧力目標値に維持されるように、前記熱浴内の圧力を制御するシミュレーション方法。
2. 前記熱浴内の圧力の制御は、前記熱浴内へ粒子を追加するか、または前記熱浴内から粒子を除去することにより行う請求項１に記載のシミュレーション方法。
3. さらに、
   前記流出入界面における温度目標値を取得し、
   粒子の状態が時間発展したときに前記熱浴内の温度が前記温度目標値に維持されるように、前記熱浴内の温度を制御する請求項１または２に記載のシミュレーション方法。
4. 流出入界面を持つ流れ場を解析領域とし、流れ場内の流体を複数の粒子の集合体として分子動力学法を用いてシミュレーションするシミュレーション装置であって、
   前記流出入界面における圧力目標値が入力される入力部と、
   前記流出入界面に、境界領域を介して熱浴を接続し、前記熱浴と前記境界領域との間、及び前記境界領域と前記解析領域との間で粒子の移動を許容した解析モデルにおいて、粒子の状態が時間発展させるとともに、粒子の状態が時間発展したときに前記境界領域内の圧力が前記圧力目標値に維持されるように、前記熱浴内の圧力を制御する処理部と
   を有するシミュレーション装置。
5. 前記熱浴内の圧力の制御は、前記熱浴内へ粒子を追加するか、または前記熱浴内から粒子を除去することにより行う請求項４に記載のシミュレーション装置。
6. 前記入力部に、さらに、前記流出入界面における温度目標値が入力され、
   前記処理部は、さらに、粒子の状態が時間発展したときに前記熱浴内の温度が前記温度目標値に維持されるように、前記熱浴内の温度を制御する請求項４または５に記載のシミュレーション装置。
7. 流出入界面を持つ流れ場を解析領域とし、流れ場内の流体を複数の粒子の集合体として分子動力学法を用いてシミュレーションする機能と、
   前記流出入界面における圧力目標値を取得する機能と、
   前記流出入界面に、境界領域を介して熱浴を接続し、前記熱浴と前記境界領域との間、及び前記境界領域と前記解析領域との間で粒子の移動を許容した解析モデルにおいて、粒子の状態が時間発展したときに前記境界領域内の圧力が前記圧力目標値に維持されるように、前記熱浴内の圧力を制御する機能と
   をコンピュータに実行させるプログラム。
8. 前記熱浴内の圧力の制御は、前記熱浴内へ粒子を追加するか、または前記熱浴内から粒子を除去することにより行う請求項７に記載のプログラム。
9. さらに、
   前記流出入界面における温度目標値を取得する機能と、
   粒子の状態が時間発展したときに前記熱浴内の温度が前記温度目標値に維持されるように、前記熱浴内の温度を制御する機能と
   をコンピュータに実行させる請求項７または８に記載のプログラム。
   

Images