JP6869621B2 - シミュレーション方法及びシミュレーション装置 - Google Patents

Description

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

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

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

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

分子動力学法を用いた従来のシミュレーション解析では、流体の流入及び流出を伴うような系を取り扱うことが困難であった。本発明の目的は、流体の流入及び流出を伴う系を、分子動力学法を用いて解析するシミュレーション方法及びシミュレーション装置を提供することである。

本発明の一観点によると、
流出入界面を持つ流れ場を解析領域とし、流れ場内の流体を複数の粒子の集合体として分子動力学法を用いてシミュレーションする方法であって、
前記解析領域の流出入界面に熱浴を接続し、前記熱浴と前記解析領域との間で粒子の移動を許容した解析モデルにおいて、時間の経過によって生じる前記熱浴内の温度と圧力との変動を補償して前記熱浴内の温度と圧力とを目標値に維持する処理を含み、
前記熱浴内の圧力を目標値に維持する処理において、
現時点における前記熱浴内の粒子の数及び状態に基づいて次状態における前記熱浴内の圧力を算出し、
算出された前記熱浴内の次状態の圧力と、圧力の目標値とを比較し、比較結果に基づいて前記熱浴内に追加または前記熱浴内から除去する粒子の個数を算出し、
算出された個数の粒子の追加または除去を前記熱浴に対して実行するシミュレーション方法が提供される。

本発明の他の観点によると、
入出力手段と、
流出入界面を持つ流れ場を解析領域とし、流れ場内の流体を複数の粒子の集合体として分子動力学法を用いて粒子の挙動を解析する処理手段と
を有し、
前記処理手段は、
前記解析領域の流出入界面の温度の目標値、圧力の目標値、及び初期条件を、前記入出力手段から取得し、
前記解析領域の流出入界面に熱浴を接続し、前記熱浴と前記解析領域との間で粒子の移動を許容した解析モデルにおいて、時間の経過によって生じる前記熱浴内の温度と圧力との変動を補償して前記熱浴内の温度と圧力とを目標値に維持しながら、前記解析領域及び前記熱浴内の粒子の挙動を、分子動力学法を用いて解析し、
前記熱浴内の圧力を目標値に維持する処理において、
現時点における前記熱浴内の粒子の数及び状態に基づいて前記熱浴内の圧力を算出し、
算出された前記熱浴内の現時点の圧力と、圧力の目標値とを比較し、比較結果に基づいて前記熱浴内に追加または除去する粒子の個数を算出し、
算出された個数の粒子の追加または除去を前記熱浴に対して実行し、
解析結果を前記入出力手段から出力させるシミュレーション装置が提供される。

解析領域の流出入界面に熱浴を接続し、熱浴内の温度と圧力とを目標値に維持することにより、流出入界面における温度と圧力とを制御して解析を行うことができる。

図１は、実施例によるシミュレーション方法で解析する対象となる解析モデルの概略図である。 図２は、実施例によるシミュレーション方法のフローチャートである。 図３Ａは、図２のステップＳＡ５の詳細なフローチャートであり、図３Ｂは、個数ｄＮが正の場合の粒子の追加の様子を示す模式図であり、図３Ｃは、個数ｄＮが負の場合の粒子の除去の様子を示す模式図である。 図４は、本変形例を説明するための第１の熱浴の模式図である。 図５は、第１の熱浴内の総粒子数の時間変化を示すグラフである。 図６Ａは、実施例によるシミュレーション方法を適用した時の熱浴内の温度と圧力との時間変化を示すグラフであり、図６Ｂは、圧力制御を行わない比較例によるシミュレーション方法を適用した時の熱浴内の温度と圧力との時間変化を示すグラフである。 図７は、実施例によるシミュレーション装置のブロック図である。

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

解析領域２０に、流出入界面２１、２２を介してそれぞれ第１の熱浴１１及び第２の熱浴１２を接続した解析モデルを定義する。第１の熱浴１１及び第２の熱浴１２の各々は例えば六面体で構成されており、流出入界面２１、２２に接する面を通って粒子３０が流出入する。このように、第１の熱浴１１及び第２の熱浴１２と解析領域２０との間で粒子３０の移動が許容されている。第１の熱浴１１及び第２の熱浴１２の他の５面には、反射境界条件を適用する。反射境界に接触した粒子３０は、接触時とは反対の面法線方向の速度成分をもって反射される。解析領域２０の壁面２３には、例えば周期境界条件または反射境界条件を適用する。第１の熱浴１１内の圧力を第２の熱浴１２内の圧力より高く設定すると、解析領域２０内に流出入界面２１から２２に向かう流れ場が形成される。

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

図２は、実施例によるシミュレーション方法のフローチャートである。まず、シミュレーションの初期条件と境界条件とを設定する（ステップＳＡ１）。初期条件は、例えば粒子数、粒子３０の質量、複数の粒子３０の位置と速度で与えられる。境界条件は、流出入界面２１、２２（図１）における温度と圧力で与えられる。境界条件で設定される温度及び圧力を、温度及び圧力の目標値ということとする。次に、初期条件に基づいて、解析領域２０、第１の熱浴１１、及び第２の熱浴１２内に複数の粒子３０を配置する（ステップＳＡ２）。

粒子３０の間の相互作用に基づいて運動方程式を解くことにより、粒子３０の次状態を算出する（ステップＳＡ３）。具体的には、１タイムステップ後における粒子３０の位置及び速度を算出する。粒子３０の間の相互作用として、例えばレナード−ジョーンズポテンシャルを用いることができる。

第１の熱浴１１及び第２の熱浴１２内の粒子３０の次状態に基づいて、第１の熱浴１１及び第２の熱浴１２内の圧力を算出する（ステップＳＡ４）。第１の熱浴１１及び第２の熱浴１２内の圧力は、例えばビリアル定理を用いて計算することができる。

ステップＳＡ４で算出された第１の熱浴１１及び第２の熱浴１２内の次状態の圧力Ｐ１、Ｐ２の、目標値からの変動を補償し、圧力Ｐ１、Ｐ２を目標値に維持する処理を行う（ステップＳＡ５）。以下、圧力を目標値に維持する処理を、圧力制御ということとする。この処理の詳細については、後に図３Ａ〜図３Ｃを参照して説明する。

次に、第１の熱浴１１及び第２の熱浴１２内の温度Ｔ１、Ｔ２を、温度の目標値に維持するための温度制御を行う（ステップＳＡ６）。温度制御には、例えば速度スケーリング法を用いることができる。

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

解析を終了するまで、ステップＳＡ３からステップＳＡ７までの処理を繰り返す（ステップＳＡ８）。解析を終了する場合には、シミュレーショ結果を出力する（ステップＳＡ９）。

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

ステップＳＡ３で算出した次状態における第１の熱浴１１内の粒子３０の総数をＮ、次状態における第１の熱浴１１内の圧力をＰｎ、圧力の目標値をＰｔで表したとき、追加すべき粒子の個数ｄＮを、以下の式で算出する。

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

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

図３Ｂは、個数ｄＮが正の場合の粒子３０の追加の様子を示す模式図である。図３Ｂにおいて、現時点の粒子３０の位置を破線で表し、次状態の粒子３０の位置を実線で表している。現時点から次状態に遷移するときに、２個の粒子３０ｃが第１の熱浴１１から流出入界面２１を通過して外側に流出している。次状態では、第１の熱浴１１内の粒子数が減少するため、圧力が低下する。このため、Ｐｔ＞Ｐｎが成立し、個数ｄＮが正になる。一例として、ｄＮ＝２の場合、図３Ｂの下図に示すように、第１の熱浴１１内に２つの粒子３０ａを追加する。

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

図３Ｃは、個数ｄＮが負の場合の粒子３０の除去の様子を示す模式図である。図３Ｃにおいて、現時点の粒子３０の位置を破線で表し、次状態の粒子３０の位置を実線で表している。現時点から次状態に遷移するときに、２個の粒子３０ｄが第２の熱浴１２に流出入界面２２を通過して流入している。次状態では、第２の熱浴１２内の粒子数が増加するため、圧力が上昇する。このため、Ｐｔ＜Ｐｎが成立し、個数ｄＮが負になる。一例として、ｄＮ＝２の場合、図３Ｃの下図に示すように、第２の熱浴１２内から２つの粒子３０ｂを除去する。

次に、上記実施例によるシミュレーション方法が有する優れた効果について説明する。
解析領域２０の流出入界面２１、２２にそれぞれ第１の熱浴１１及び第２の熱浴１２を接続して、第１の熱浴１１及び第２の熱浴１２の温度制御及び圧力制御を行うことにより、解析領域２０の流出入界面２１、２２における温度及び圧力を制御することができる。

また、上記実施例では、粒子数の増減に起因する圧力の変動に応じて、圧力が目標値に近づくように粒子の追加または除去を行っている。このため、解析の開始から一定の時間が経過すると、第１の熱浴１１及び第２の熱浴１２内の圧力及び粒子数がほぼ目標値に維持された平衡状態に達する。これにより、圧力を一定の目標値に維持することができる。第１の熱浴１１及び第２の熱浴１２内の圧力制御を行うことにより、解析領域２０の流出入界面２１、２２における圧力を目標値に維持することができる。

次に、図４を参照して、上記実施例の変形例によるシミュレーション方法について説明する。上記実施例では、第１の熱浴１１内の圧力が低下した時、第１の熱浴１１内の任意の位置に新たな粒子３０ａ（図３Ｂ）を追加し、追加した粒子３０ａと既存の粒子３０との間に作用する斥力が大きくなりすぎる場合には、粒子を再配置した。本変形例では、追加した粒子３０ａと既存の粒子３０との間に作用する斥力が大きくなりすぎる事態の発生確率を低減させる。

図４は、本変形例を説明するための第１の熱浴１１の模式図である。本変形例では、第１の熱浴１１内を体積の等しい複数のセル４０に区分する。なお、セル４０は、仮想的に区分されているのであり、セル４０の境界は粒子３０の移動に何ら影響を及ぼさない。図４では、第１の熱浴１１のある断面におけるセル４０を示している。実際には、第１の熱浴１１は、縦横高さ方向に三次元的に区分される。ステップＳＢ２（図３Ａ）の前に、セル４０ごとに次状態における粒子３０の数ｎを算出する。図４では、一例として、右上のセル４０の粒子数ｎが２、中央列の最上段のセル４０の粒子数ｎが３、その他のセル４０の粒子数が４である例を示している。

ステップＳＢ２において、相対的に粒子数が少ないセル４０に優先的に新たな粒子３０を追加する。例えば、まず、最も粒子数の少ないセル４０（図４において右上のセル４０）に、新たな粒子３０ｅを追加する。粒子が追加されたセル４０の粒子数は、１個分増加する。次に、その時点で最も粒子数が少ないセル４０（図４において右上または中央列の最上段のセル４０）に、新たな粒子３０ｆを追加する。追加すべき個数ｄＮの粒子が追加されるまで、この処理を繰り返す。

本変形例では、第１の熱浴１１内の、粒子の密度が相対的に低い箇所に新たな粒子３０ａを追加するため、追加した粒子３０ａと既存の粒子３０との間に作用する斥力が大きくなりすぎる事態の発生確率を低減させることができる。

第２の熱浴１２内の粒子３０を除去する場合には、粒子数が相対的に多いセル４０から除去すべき粒子３０ｂ（図３Ｃ）を抽出するとよい。

上記実施例では、第１の熱浴１１を流れ場の上流側の流出入界面２１に接続し、第２の熱浴１２を流れ場の下流側の流出入界面２２に接続した。第１の熱浴１１及び第２の熱浴１２のいずれが上流側となるかは、流出入界面２１、２２の境界条件によって決定されるのであり、シミュレーションにおいて第１の熱浴１１と第２の熱浴１２とを区別して取り扱う必要はない。また、解析のタイムステップによっては、短時間のゆらぎによって、上流側に接続された第１の熱浴１１から粒子を除去し、下流側に接続された第２の熱浴１２に新たな粒子を追加する場合も生じ得る。

また、上記実施例では、タイムステップごとに圧力制御を行ったが、必ずしもすべてのタイムステップにおいて圧力制御を行う必要は無い。例えば、所定の複数のタイムステップごとに圧力制御を行ってもよい。または、第１の熱浴１１及び第２の熱浴１２内の圧力と、圧力の目標値との差が所定の許容上限値を超えた時に圧力制御を行うようにしてもよい。

次に、図５Ａ〜図６Ｂを参照して、第１の熱浴１１についてシミュレーションを行った結果について説明する。

図５は、第１の熱浴１１内の総粒子数の時間変化を示すグラフである。横軸はタイムステップを表し、縦軸は総粒子数を表す。図５のグラフにおいて、実線は圧力制御を行う実施例によるシミュレーション方法を適用した結果を示し、破線は圧力制御を行わない比較例によるシミュレーション方法を適用した結果を示す。

圧力制御を行わない比較例においては、第１の熱浴１１から粒子が流出するため、時間の経過とともに総粒子数が減少する。圧力制御を行う実施例においては、タイムステップ数が約３００に達するまでは、総粒子数が時間の経過に伴って減少するが、その後は総粒子数がほぼ一定に維持される。最初に総粒子数が減少するのは、熱浴内の圧力の初期値を、圧力の目標値より高く設定したためである。タイムステップが約３００まで進むと、熱浴内の圧力が目標値にほぼ等しくなる。実施例では、圧力を目標値に維持するように、新規の粒子の追加が行われるため、総粒子数は、ほぼ一定になる。

図６Ａは、実施例によるシミュレーション方法を適用した場合の熱浴内の温度と圧力との時間変化を示すグラフである。図６Ｂは、圧力制御を行わない比較例によるシミュレーション方法を適用した場合の熱浴内の温度と圧力との時間変化を示すグラフである。図６Ａ及び図６Ｂの横軸はタイムステップを表し、縦軸は温度及び圧力を任意単位で表す。グラフ中の太い実線及び細い実線が、それぞれ圧力及び温度の時間変化を示す。なお、実施例及び比較例のいずれにおいても、温度制御のアルゴリズムにより温度が目標値に維持されるように制御されている。

実施例（図６Ａ）では、タイムステップ数が５００を超えた時点から圧力がほぼ一定に維持されている。これに対し、比較例（図６Ｂ）では、タイムステップが５００を超えても時間の経過とともに圧力が減少している。これは、熱浴内の粒子が流出入界面を通って外部に流出するためである。実施例では、熱浴内に新たに粒子が追加されるため、圧力が目標値に維持される。

実施例及び比較例のいずれの場合でも温度制御を行っているため、時間が経過するに従って温度が目標値に近づき、その後温度はほぼ目標値に維持される。

次に、図７を参照して、上記実施例によるシミュレーション方法の各工程を実行するシミュレーション装置について説明する。

図７は、実施例によるシミュレーション装置のブロック図である。シミュレーション装置は一般的なコンピュータで構成され、中央処理ユニット（ＣＰＵ、処理手段）５０、メモリ５１、入出力装置（入出力手段）５２、及び補助記憶装置５３を含む。これらは、バス５４で相互に接続されている。

入出力装置５２は、キーボード、マウス等のポインティングデバイス、ディスプレイ、リムーバブルメディアのリーダライタ、通信装置等を含む。ディスプレイは、ユーザの操作に必要な各種ウィンドウやデータ等を表示する。通信装置は、外部機器とデータ通信を行う。ユーザは、入出力装置を操作することにより、コンピュータに対する指示、各処理に必要なデータの入力等を行う。また、処理結果が、入出力装置に出力される。

補助記憶装置５３は、ハードディスク等で構成され、コンピュータのＯＳ、シミュレーションプログラム、シミュレーションに必要なデータや処理結果等を蓄積する。

メモリ５１は、ＲＯＭ（Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）やＲＡＭ（Ｒａｎｄｏｍ Ａｃｃｅｓｓ Ｍｅｍｏｒｙ）等で構成される。メモリ５１は、ＣＰＵ５０により補助記憶装置５３から読み出されたシミュレーションプログラム等を格納する。

ＣＰＵ５０は、コンピュータのＯＳ、メモリ５１に格納されているシミュレーションプログラム等に基づいて、各種演算や入出力装置５２及び補助記憶装置５３の制御を行う。

次に、シミュレーション装置の動作について説明する。ＣＰＵ５０は、ステップＳＡ１（図２）の処理を実行するために、初期条件及び境界条件を入出力装置５２から取得する。例えば、ユーザが入出力装置５２を操作することにより、初期条件及び境界条件を入力する。または、ＣＰＵ５０は、初期条件及び境界条件をリムーバブルメディアから読み込む。その後、ＣＰＵ５０は、ステップＳＡ２を実行し、ステップＳＡ３からステップＳＡ８までの処理を繰り返す。ステップＳＡ９では、ＣＰＵ５０は、例えば解析結果をディスプレイに表示させる。

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

１１ 第１の熱浴
１２ 第２の熱浴
２０ 解析領域
２１、２２ 流出入界面
２３ 壁面
３０ 粒子
３０ａ 新たに追加された粒子
３０ｂ 除去される粒子
４０ 第１の熱浴を区分したセル
５０ 中央処理ユニット（ＣＰＵ）
５１ メモリ
５２ 入出力装置
５３ 補助記憶装置
５４ バス

Claims (4)

1. 流出入界面を持つ流れ場を解析領域とし、流れ場内の流体を複数の粒子の集合体として分子動力学法を用いてシミュレーションする方法であって、
   前記解析領域の流出入界面に熱浴を接続し、前記熱浴と前記解析領域との間で粒子の移動を許容した解析モデルにおいて、時間の経過によって生じる前記熱浴内の温度と圧力との変動を補償して前記熱浴内の温度と圧力とを目標値に維持する処理を含み、
   前記熱浴内の圧力を目標値に維持する処理において、
   現時点における前記熱浴内の粒子の数及び状態に基づいて次状態における前記熱浴内の圧力を算出し、
   算出された前記熱浴内の次状態の圧力と、圧力の目標値とを比較し、比較結果に基づいて前記熱浴内に追加または前記熱浴内から除去する粒子の個数を算出し、
   算出された個数の粒子の追加または除去を前記熱浴に対して実行するシミュレーション方法。
2. 粒子の追加または除去を前記熱浴に対して実行する処理において、
   前記熱浴の領域内を区分して得られる複数のセルごとに前記セル内の粒子の個数を算出し、相対的に粒子数の少ない前記セルに優先的に粒子を追加し、相対的に粒子数の多い前記セルから優先的に粒子を除去する請求項１に記載のシミュレーション方法。
3. 入出力手段と、
   流出入界面を持つ流れ場を解析領域とし、流れ場内の流体を複数の粒子の集合体として分子動力学法を用いて粒子の挙動を解析する処理手段と
   を有し、
   前記処理手段は、
   前記解析領域の流出入界面の温度の目標値、圧力の目標値、及び初期条件を、前記入出力手段から取得し、
   前記解析領域の流出入界面に熱浴を接続し、前記熱浴と前記解析領域との間で粒子の移動を許容した解析モデルにおいて、時間の経過によって生じる前記熱浴内の温度と圧力との変動を補償して前記熱浴内の温度と圧力とを目標値に維持しながら、前記解析領域及び前記熱浴内の粒子の挙動を、分子動力学法を用いて解析し、
   前記熱浴内の圧力を目標値に維持する処理において、
   現時点における前記熱浴内の粒子の数及び状態に基づいて前記熱浴内の圧力を算出し、
   算出された前記熱浴内の現時点の圧力と、圧力の目標値とを比較し、比較結果に基づいて前記熱浴内に追加または除去する粒子の個数を算出し、
   算出された個数の粒子の追加または除去を前記熱浴に対して実行し、
   解析結果を前記入出力手段から出力させるシミュレーション装置。
4. 前記処理手段は、粒子の追加または除去を前記熱浴に対して実行する処理において、
   前記熱浴の領域内を区分して得られる複数のセルごとに前記セル内の粒子の個数を算出し、相対的に粒子数の少ない前記セルに優先的に粒子を追加し、相対的に粒子数の多い前記セルから優先的に粒子を除去する請求項３に記載のシミュレーション装置。
   

Images