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

Description

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

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

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

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

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

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

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

解析領域の流出入界面に熱浴を接続し、熱浴内の圧力制御を行うことにより、流出入界面における圧力を目標値に維持した状態でシミュレーションを行うことができる。

図１は、実施例によるシミュレーション装置のブロック図である。 図２Ａは、実施例によるシミュレーション方法で解析する対象となる解析モデルの一例を示す概略図であり、図２Ｂは、流出入界面の面内位置と、圧力目標値との関係の一例を示すグラフである。 図３は、実施例によるシミュレーション装置の処理部が実行する処理手順のフローチャートである。 図４は、流出入界面における圧力制御の手順を示すフローチャートである。 図５Ａは、個数ｄＮが正の場合の粒子の追加の様子を示す第１熱浴セルの模式図であり、図５Ｂは、個数ｄＮが負の場合の粒子の除去の様子を示す第２熱浴の模式図である。 図６は、実際に行ったシミュレーションの解析モデルの模式図である。 図７Ａに、流出入界面における圧力目標値と、定常状態に達したときの第１熱浴内の圧力の計算値を示すグラフであり、図７Ｂは、定常状態に達したときの粒子の位置を示す図である。 図８は、他の実施例によるシミュレーション方法で解析する対象となる解析モデルの斜視図である。 図９は、圧力制御（図３のステップＳ４）を行う手順のフローチャートである。

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

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

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

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

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

流出入界面２１が複数の区画２１ａに区分されている。流出入界面２１は、例えば格子状に区分される。第１熱浴３０が、複数の区画２１ａに対応して複数の第１熱浴セル３０ａに区分されている。複数の第１熱浴セル３０ａは、それぞれ対応する区画２１ａを介して解析領域２０に接続される。隣接する２つの区画２１ａに対応する２つの第１熱浴セル３０ａは相互に連続しており、両者の界面は、２つの第１熱浴セル３０ａの間で粒子５０の移動を許容する。

図２Ｂは、流出入界面２１の面内位置と、圧力目標値Ｐ１ｏとの関係の一例を示すグラフである。横軸は、流出入界面２１の面内の一次元方向の位置を表し、縦軸は圧力を表す。流出入界面２１の面内に関して分布を持つ圧力目標値Ｐ１ｏが、シミュレーション条件の１つとして与えられる。圧力目標値Ｐ１ｏの分布を反映するように、区画２１ａごとに圧力目標値Ｐ１ａｏが決定される。例えば、圧力目標値Ｐ１ｏを区画２１ａの寸法で離散化し、区画２１ａにおける圧力目標値Ｐ１ａｏとして、離散化した値を設定する。

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

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

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

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

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

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

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

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

流出入界面２１における温度目標値Ｔ１ｏが面内で一定ではなく、分布を持っている場合には、温度制御を第１熱浴セル３０ａごとに行う。

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

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

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

次に、現時点のタイムステップにおいて圧力制御を実行するか否かを判定する（ステップＳ４２）。例えば、圧力制御は、数百タイムステップごとに１回行う。圧力制御を実行しない場合は、圧力制御の処理を終了して、図３に示したフローチャートに戻る。圧力制御を実行する場合は、第１熱浴セル３０ａ内の圧力Ｐ１ｈａの平均値と、対応する区画２１ａにおける圧力目標値Ｐ１ａｏとに基づいて、第１熱浴セル３０ａ内の圧力が圧力目標値Ｐ１ａｏに維持されるように、第１熱浴セル３０ａに対して追加または除去する粒子数を算出する（ステップＳ４３）。第１熱浴セル３０ａ内の圧力Ｐ１ｈａの平均値は、ステップＳ４１で求められた過去の複数のタイムステップにおける圧力Ｐ１ｈａを平均することにより求められる。

以下、追加または除去する粒子数の計算方法について説明する。
ステップＳ３で算出した最新の状態における第１熱浴セル３０ａ内の粒子５０の総数をＮとしたとき、追加すべき粒子の個数ｄＮを、以下の式で算出する。

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

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

次に、流出入界面２２（図２）における圧力制御（図３のステップＳ４）について説明する。流出入界面２２においては、第２熱浴４０内の圧力Ｐ２ｈの平均値と、圧力目標値Ｐ２ｏを用いて第２熱浴４０内の圧力を制御する。すなわち、第２熱浴４０に対しても、第１熱浴セル３０ａと同様に粒子５０の追加または除去を行う。

第２熱浴４０については、第２熱浴４０内の粒子５０の総数をＮとしたとき、追加すべき粒子の個数ｄＮは、以下の式で算出される。

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

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

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

次に、上記実施例の優れた効果について説明する。
上記実施例では、第１熱浴セル３０ａ内の圧力Ｐ１ｈａが、流出入界面２１の対応する区画２１ａにおける圧力目標値Ｐ１ａｏに維持されるため、流出入界面２１における圧力境界条件が満足される。同様に、もう一方の流出入界面２２における圧力境界条件が満足される。このように、圧力境界条件が満たされる状態で、シミュレーションを行うことができる。

さらに、上記実施例では、流出入界面２１を複数の区画２１ａに区分しているため、流出入界面２１の面内において圧力目標値が一様でない場合でも、圧力境界条件が満たされる状態でシミュレーションを行うことができる。

次に、上記実施例の変形例について説明する。上記実施例では、流出入界面２１を複数の区画２１ａに区分し、もう一方の流出入界面２２については、圧力目標値Ｐ２ｏが面内で一様であると仮定している。流出入界面２２における圧力目標値Ｐ２ｏが面内で一様ではない場合に、流出入界面２２も複数の区画に区分し、区画ごとに圧力目標値を設定するとよい。

また、上記実施例では、解析領域２０の壁面２３（図２Ａ）を反射境界または周期境界としてシミュレーションを行っているが、壁面２３に圧力境界条件を設定してもよい。この場合には、壁面２３に熱浴を接続するとよい。圧力が面内で一様でない場合には、壁面２３を複数の区画に区分するとよい。

上記実施例では、温度制御について詳細な説明は省略したが、流出入界面２１の面内方向に温度分布が発生している温度境界条件を適用する場合には、流出入界面２１の区画２１ａ及び第１熱浴セル３０ａごとに温度制御を行うとよい。

次に、図６～図７Ｂを参照して、上記実施例による効果を確認するために行ったシミュレーションについて説明する。

図６は、実際に行ったシミュレーションの解析モデルの模式図である。解析領域２０の流出入界面２１に第１熱浴３０を接続した。流出入界面２１とは反対側においては、解析領域２０が開放されている。流出入界面２１を一方向に１５等分して１５個の区画２１ａを定義した。第１熱浴３０も、区画２１ａに対応して１５等分し、１５個の第１熱浴セル３０ａを定義した。

図７Ａに、流出入界面２１の区画２１ａにおける圧力目標値Ｐ１ａｏと、定常状態に達したときの第１熱浴３０内の圧力の計算値を示すグラフである。図７Ａの横軸は、流出入界面２１の面内位置を表す。横軸の両端が、それぞれ流出入界面２１の両端に相当する。縦軸は、圧力を無次元量として表す。図７Ａに示した実線は、シミュレーション条件として与えた圧力目標値Ｐ１ａｏを示し、丸記号は、シミュレーションにより求められた第１熱浴セル３０ａ内の圧力Ｐ１ｈａを示す。圧力目標値Ｐ１ａｏに、中央部分が両端部分より相対的に高くなるような分布を与えた。第１熱浴セル３０ａ内の圧力Ｐ１ｈａが圧力目標値Ｐ１ａｏを反映した分布を有していることがわかる。

図７Ｂは、定常状態に達したときの粒子の位置を示す図である。複数の第１熱浴セル３０ａ内の粒子密度が、中央の第１熱浴セル３０ａから両端の第１熱浴セル３０ａに向かって徐々に低くなっていることが確認される。区画２１ａにおける圧力目標値Ｐ１ａｏの分布が、実際の粒子密度の分布に反映されていることがわかる。

図６～図７Ｂに示したシミュレーション結果から、複数の粒子の状態を時間発展させたときに、シミュレーション条件として与えた圧力境界条件が維持されていることが確認された。

次に、図８及び図９を参照して、他の実施例について説明する。以下、図１～図５Ｂに示した実施例と共通の構成については説明を省略する。

図８は、本実施例によるシミュレーション方法で解析する対象となる解析モデルの斜視図である。図２Ａに示した実施例では、第１熱浴３０が解析領域２０に直接接しているが、本実施例では、第１熱浴３０が境界領域３１を介して解析領域２０に接続されている。境界領域３１も、流出入界面２１の複数の区画２１ａに対応して複数の境界領域セル３１ａに区分されている。

本実施例においては、境界領域セル３１ａ内の圧力が、対応する区画２１ａの圧力目標値Ｐ１ａｏに維持されるように、第１熱浴セル３０ａ内の圧力Ｐ１ｈａを制御する。境界領域セル３１ａに対しては、粒子の追加や除去を行わない。

図９は、圧力制御（図３のステップＳ４）を行う手順のフローチャートである。
まず、複数の境界領域セル３１ａ内の圧力を計算する（ステップＳ４５）。境界領域セル３１ａ内の圧力は、例えばビリアル定理を用いて計算することができる。

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

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

式（３）において、左辺のＰ１ｈａｏは更新後の圧力目標値であり、右辺のＰ１ｈａｏは現時点（更新前）の圧力目標値である。すなわち、流出入界面２１の区画２１ａにおける圧力目標値Ｐ１ａｏと境界領域セル３１ａの圧力Ｐ１ｂａの平均値との比、及び、第１熱浴セル３０ａの現時点の圧力目標値Ｐ１ｈａｏに基づいて、圧力目標値Ｐ１ｈａｏを更新する。

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

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

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

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

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

１０ 入力部
１１ 処理部
１２ 出力部
１３ 記憶部
２０ 解析領域
２１ 流出入界面
２１ａ 区画
２２ 流出入界面
２３ 壁面
３０ 第１熱浴
３０ａ 第１熱浴セル
３１ 境界領域
３１ａ 境界領域セル
４０ 第２熱浴
５０、５０ａ、５０ｂ、５０ｃ、５０ｄ 粒子

Claims (9)

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

Images