JP5930987B2 - 解析装置およびコンピュータプログラム - Google Patents

Description

本発明は、粒子系を解析する解析装置に関する。

従来、古典力学や量子力学等を基に計算機を用いて物質科学全般の現象を探るための方法として、分子動力学法（Molecular Dynamics Method、以下ＭＤ法と称す）に基づくシミュレーションが知られている。ＭＤ法に基づくシミュレーションでは、粒子の運動を物理的により厳密に取り扱えるため相転移等も好適に扱える反面、マクロなサイズを扱うと計算量も飛躍的に増大するので実用的にはより少ない数の粒子しか扱えない。したがって従来ＭＤ法は主に材料の物性の予測等、ミクロなスケールや解析対象の形状があまり関係しない用途に使用されることが多かった。

最近、このＭＤ法をマクロスケールの系を扱えるように発展させた繰り込み群分子動力学法（Renormalized Molecular Dynamics、以下ＲＭＤ法と称す）が提案された（例えば、特許文献１参照）。

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

気液二相流や固気二相流などの混相流の界面に特徴的な現象をＲＭＤ法により解析する場合において、計算時間の短縮のために繰り込み階数を増すと、密度が低い方の相の粒子数が極端に少なくなり連続体として振る舞わなくなる場合がある。すなわち、その相が希薄流れとなる場合がある。この場合、解析結果は実際とは乖離する。

本発明はこうした課題に鑑みてなされたものであり、その目的は、混相流の界面に特徴的な現象をより正確に扱える解析技術の提供にある。

上記課題を解決するために、本発明のある態様の解析装置は、繰り込み群分子動力学法を使用して混相流を解析する解析装置であって、混相流を記述する複数の粒子を含む粒子系を取得する粒子系取得部と、粒子間の距離に基づき粒子に働く力を演算する力演算部と、離散化された粒子の運動方程式に力演算部によって演算された力を適用することによって粒子の位置および速度のうちの少なくともひとつを演算する粒子状態演算部と、を備える。混相流は、少なくとも、第１相と当該第１相より粒子密度が低い第２相とが混在した流れであり、第１相の粒子には原子間ポテンシャルによって決定される粒子径が設定され、第２相の粒子には第２相の原子間ポテンシャルによって決定される粒子径よりも大きい粒子径が設定されている。

この態様によると、粒子密度が低い第２相が連続体として振る舞う。

なお、以上の構成要素の任意の組み合わせや、本発明の構成要素や表現を装置、方法、システム、コンピュータプログラム、コンピュータプログラムを格納した記録媒体などの間で相互に置換したものもまた、本発明の態様として有効である。

本発明によれば、混相流の界面に特徴的な現象をより正確に扱うことができる。

本実施の形態の解析対象を示す模式図である。 本実施の形態に係る解析装置の機能および構成を示すブロック図である。 図２の解析装置における一連の処理の一例を示すフローチャートである。 図４（ａ）〜（ｔ）は、本実施の形態に係る手法を使用した演算結果を示す図である。 図５（ａ）〜（ｔ）は、従来の手法を使用した演算結果を示す図である。 実施の形態の変形例の解析対象を示す模式図である。

以下、本発明を好適な実施の形態をもとに図面を参照しながら説明する。各図面に示される同一または同等の構成要素、部材、処理には、同一の符号を付するものとし、適宜重複した説明は省略する。

以下、本実施の形態の原理を説明する。
図１は本実施の形態の解析対象の粒子系を示す模式図である。図１では、液相２が仮想の壁３にせき止められており、液相２の周りはその蒸気である気相４で満たされている。本実施の形態では、図１において、壁３を取り外した場合に起こる現象をＲＭＤ法により解析することを前提に説明する。つまり、気液二相流の界面に特徴的な現象をＲＭＤ法により解析することを前提に説明する。

ＲＭＤ法により解析する場合、ユーザはまず、解析装置や他の演算装置を使用して解析対象を３次元の仮想空間内に生成する。ユーザは、仮想空間内にＭＤ法における粒子すなわち現実世界の原子または分子に対応する粒子を複数配置する。一例では、配置する粒子の数はアボガドロ数程度となる。

以下では粒子は全て同質または同等なものとして設定され、かつ、ポテンシャルエネルギ関数は２体のポテンシャルであって粒子によらずに同じ形を有するものとして設定される場合について説明する。しかしながら、他の場合にも本実施の形態に係る技術的思想を適用できることは、本明細書に触れた当業者には明らかである。

続いてユーザは、配置された複数の粒子からなる系に、ＲＭＤ法における変換則を適用する。一例では、数十回の繰り込みにより、粒子数は数万個程度まで低減される。

ここで、繰り込み階数を増すと、液相２に比べて密度が低い気相４の粒子数は極端に減少し、ついには気相４は連続体として振る舞わなくなる。このことは、繰り込みに際し平均自由行程λが不変であるため、繰り込まれた系のクヌーセン数Ｋｎ′が増大することからも分かる。クヌーセン数Ｋｎ′は以下の式（１）により与えられる。
ここで、Ｌ′は系の代表長さ、αはくりこみ因子、ｎは繰り込み回数である。

気相４が連続体として振る舞うようにするためには、繰り込み際して気相４の粒子（以下、「気体粒子」ともよぶ）の平均自由行程λも同時に変換する、すなわち平均自由行程λを小さくしてクヌーセン数Ｋｎ′を小さくしてすればよい。平均自由行程λは以下の式（２）により与えられる。
ここで、ｎは濃度、ｄは粒子径である。

濃度ｎは、温度と圧力が指定されれば状態方程式（例えば、理想気体ならＰ＝ｎｋ_ＢＴ）から決まる。このため、平均自由行程λを小さくしてクヌーセン数Ｋｎ′を小さくするには粒子径を増大させればよい。

したがって、ユーザは、繰り込まれた系に含まれる粒子のうち、気体粒子の粒子径を、原子の直径（ポテンシャルで決まる径）よりも大きい粒子径に置き換えることにより、気相４を連続体として取り扱うことができ、気相４と液相２との界面に特徴的な現象をより正確に解析することができる。なお、液相２の粒子（以下、「液体粒子」ともよぶ）の粒子径は、原子の直径のままでよい。

次に、気体粒子の粒子径を決定する手順を示す。
気相４が連続体として振る舞うためには、気相４の粒子は以下の式（３）を満たす必要がある。
この式（３）と式（２）とから粒子径の下限が決まる。

また、気体相の各粒子が重なり合わないようにする必要があるため、気相４の粒子は以下の式（４）を満たす必要がある。
この式（４）から粒子径の上限が決まる。ユーザは、少なくともこれらの範囲に収まるように気体粒子の粒子径を決めればよい。

また、ヤング率Ｙ、粘度ηには、ポテンシャルパラメータε、
と以下の式（５）、（６）で与えられる関係がある。
ここで
はレナード・ジョーンズ（Lennard-Jones）型のポテンシャルエネルギ関数の２階微分である。
ここでｍは質量、ｋ_Ｂはボルツマン定数、Ｔは絶対温度である。

したがって、ヤング率Ｙ、粘度ηの値が変わらないように質量ｍやポテンシャルパラメータε、σ等を調整することで、気相４の粒子径を大きくしたことによるヤング率Ｙまたは／および粘度ηの変動を吸収できる。

次に、気体粒子の粒子系を大きくすることによる計算効率の低下を抑える手法について説明する。
ｄｔを液相２のＲＭＤ計算を実行する時間刻みとすると、以下の式（７）を満たすと計算効率がよい。
ここで、ｍ_liquidは液体粒子の密度、ｄ_liquidは液体粒子の粒子径、ε_liquid、σ_liquidは液体粒子のポテンシャルパラメータである。

レナード・ジョーンズ型のポテンシャルエネルギ関数の場合、式（３）および式（４）を満足するものとして、例えば
と変換される。ここで、３．５σ_liquidは液相２の粒子間の相互作用に対するカットオフ半径である。

したがって、式（７）および式（８）より式（９）の関係が成立する。

これを満たすように気相４の粒子密度を設定すれば、時間刻みを大きくでき、計算効率が上がる。一方、時間刻みを液体のｄｔより小さく取ることを厭わなければ、すなわち式（７）を考慮しなければ、任意の密度を設定することができる。つまり、計算効率を優先すると、密度に制限ができるため式（６）において粘度ηを所望の値に合わせらない場合がある。一方、粘度ηを所望の値に合わせることを優先すると、時間刻みが小さくなり計算効率が落ちる。両者にはトレードオフの関係がある。

次に、気体粒子の粒子系を大きくすることによる計算効率の低下を抑える別の手法について説明する。
気体粒子の粒子径を大きくすると、その粒子に相互作用が働く距離も大きくなる。そのため、粒子系を大きくすると、例えば、相互作用が働く粒子として液相２の粒子を大量に含む場合が生じうる。これにより、計算対象が増加し、計算時間が長くなる。そこで、式（１０）のように気相４の粒子に及ぶ引力をゼロとする。つまり、カットオフ半径を小さくする。なお、引力をカットするので気相４は理想気体となる。なお、これに加えて気相４の粒子径を、液相２の粒子間の相互作用に対するカットオフ半径と同じ３．５σ_liquidにすれば、従来法と同程度の計算回数を維持できる。

まとめると、解析装置１００を使用するユーザは、以下を考慮して繰り込まれた系に含まれる粒子に関する各パラメータを設定すればよい。
（i）気体粒子の粒子径を、式（３）および式（４）を満足する粒子径に置き換えることにより、平均自由行程λが小さくなりクヌーセン数Ｋｎ′が０．１６以下となる。すなわち、気相４は連続体として振る舞うことができる。
（ii）式（５）および式（６）を満たすように気体粒子の質量ｍやポテンシャルパラメータε、σ等を調整することで、気体粒子の粒子径を大きくことによるヤング率Ｙまたは／および粘度ηの変動を吸収しうる。
（iii）式（９）を満たすように気体粒子の密度を設定すれば、時間刻みを大きくでき、計算効率が上がる。
（iiii）式（１０）を満たすように、カットオフ半径を設定することで、計算対象が減り、計算時間を短縮できる。

なお、式（１）〜式（１０）により繰り込まれた系に含まれる粒子に関する各パラメータを設定する際に必要な演算は、後述の数値演算部１２０によって実施されてもよい。

次に、以上を考慮して気体粒子に関する各パラメータを設定する一例を示す。
式（３）および式（４）とを満たし、かつ、従来法と同程度の計算回数を維持するため、気体粒子の粒子径は
とした。
ポテンシャルパラメータε′は、式（５）および式（６）を考慮せずに、すなわちヤング率Ｙおよび粘度ηは調整せずに、
とする。

これより、計算効率を優先すると、気体粒子の質量は式（９）から
となる。

このとき気体相の密度は、
になる。

Ｌ′＝５０［ｎｍ］とすると、クヌーセン数は
となり気相４は連続体として振る舞う。

次に、以上を考慮して粒子に関する各パラメータが設定された粒子系を解析する場合について説明する。
図２は、解析装置１００の機能および構成を示すブロック図である。ここに示す各ブロックは、ハードウエア的には、コンピュータのＣＰＵ（central processing unit）をはじめとする素子や機械装置で実現でき、ソフトウエア的にはコンピュータプログラム等によって実現されるが、ここでは、それらの連携によって実現される機能ブロックを描いている。したがって、これらの機能ブロックはハードウエア、ソフトウエアの組合せによっていろいろなかたちで実現できることは、本明細書に触れた当業者には理解されるところである。

解析装置１００は入力装置１０２および出力装置１０４と接続される。入力装置１０２は、解析装置１００で実行される処理に関係するユーザの入力を受けるためのキーボード、マウスなどであってもよい。入力装置１０２は、インターネットなどのネットワークやＣＤ、ＤＶＤなどの記録媒体から入力を受けるよう構成されていてもよい。出力装置１０４は、ディスプレイなどの表示機器やプリンタなどの印刷機器であってもよい。

解析装置１００は、粒子系取得部１１０と、数値演算部１２０と、表示制御部１３０と、粒子データ保持部１５０と、を備える。

粒子系取得部１１０は、入力装置１０２を介してユーザから取得する入力情報に基づき、１、２または３次元の仮想空間内に定義されるＮ（Ｎは自然数）個の粒子からなる粒子系のデータを取得する。粒子系はＲＭＤ法を使用して繰り込まれた粒子系である。

粒子系取得部１１０は、入力情報に基づき仮想空間内にＮ個の粒子を配置し、配置されたそれぞれの粒子に速度を付与する。また、粒子系取得部１１０は、入力情報または粒子データ保持部１５０に記憶されている情報から粒子の質量などの以降の演算に必要なパラメータを取得する。粒子系取得部１１０は、配置された粒子の位置と、その粒子の速度と、その粒子の質量と、を対応付けて粒子データ保持部１５０に登録する。

数値演算部１２０は、粒子データ保持部１５０によって保持されるデータが表す粒子系の各粒子の運動を支配する支配方程式を数値的に演算する。特に数値演算部１２０は、離散化された粒子の運動方程式にしたがった繰り返し演算を行う。
数値演算部１２０は、力演算部１２２と、粒子状態演算部１２４と、状態更新部１２６と、終了条件判定部１２８と、を含む。

力演算部１２２は粒子データ保持部１５０によって保持される粒子系のデータを参照し、各粒子系の各粒子について、粒子間の距離に基づきその粒子に働く力を演算する。力演算部１２２は、第１粒子系の演算対象の粒子について、その演算対象の粒子との距離が所定のカットオフ距離よりも小さな粒子（以下、近接粒子と称す）を決定する。

力演算部１２２は、各近接粒子について、その近接粒子と演算対象の粒子との間のポテンシャルエネルギ関数およびその近接粒子と演算対象の粒子との距離に基づいて、その近接粒子が演算対象の粒子に及ぼす力を演算する。特に力演算部１２２は、その近接粒子と演算対象の粒子との距離の値におけるポテンシャルエネルギ関数のグラジエント（Gradient）の値から力を算出する。力演算部１２２は、近接粒子が演算対象の粒子に及ぼす力を全ての近接粒子について足し合わせることによって、演算対象の粒子に働く力を算出する。

粒子状態演算部１２４は粒子データ保持部１５０に保持される粒子系のデータを参照し、各粒子系の各粒子について、離散化された粒子の運動方程式に力演算部１２２によって演算された力を適用することによって粒子の位置および速度のうちの少なくともひとつを演算する。本実施の形態では、粒子状態演算部１２４は粒子の位置および速度の両方を演算する。

粒子状態演算部１２４は、力演算部１２２によって演算された力を含む離散化された粒子の運動方程式から粒子の速度を演算する。粒子状態演算部１２４は、第１粒子系の粒子について、蛙跳び法やオイラー法などの所定の数値解析の手法に基づき所定の微小な時間刻みΔｔを使用して離散化された粒子の運動方程式に、力演算部１２２によって演算された力を代入することによって、粒子の速度を演算する。この演算には以前の繰り返し演算のサイクルで演算された粒子の速度が使用される。

粒子状態演算部１２４は、演算された粒子の速度に基づいて粒子の位置を算出する。粒子状態演算部１２４は、第１粒子系の粒子について、所定の数値解析の手法に基づき時間刻みΔｔを使用して離散化された粒子の位置と速度の関係式に、演算された粒子の速度を適用することによって、粒子の位置を演算する。この演算には以前の繰り返し演算のサイクルで演算された粒子の位置が使用される。

状態更新部１２６は、粒子データ保持部１５０に保持される各粒子系の各粒子の位置および速度のそれぞれを、粒子状態演算部１２４によって演算された位置および速度で更新する。

終了条件判定部１２８は、数値演算部１２０における繰り返し演算を終了すべきか否かを判定する。繰り返し演算を終了すべき終了条件は、例えば繰り返し演算が所定の回数行われたことや、外部から終了の指示を受け付けたことや、粒子系が定常状態に達したことである。終了条件判定部１２８は、終了条件が満たされる場合、数値演算部１２０における繰り返し演算を終了させる。終了条件判定部１２８は、終了条件が満たされない場合、処理を力演算部１２２に戻す。すると力演算部１２２は、状態更新部１２６によって更新された粒子の位置で再び接触粒子対を特定する。

表示制御部１３０は、粒子データ保持部１５０に保持されるデータが表す各粒子系の各粒子の位置、速度に基づき、出力装置１０４に粒子系の時間発展の様子やある時刻における状態を表示させる。この表示は、静止画または動画の形式で行われてもよい。

上述の実施の形態において、保持部の例は、ハードディスクやメモリである。また、本明細書の記載に基づき、各部を、図示しないＣＰＵや、インストールされたアプリケーションプログラムのモジュールや、システムプログラムのモジュールや、ハードディスクから読み出したデータの内容を一時的に記憶するメモリなどにより実現できることは本明細書に触れた当業者には理解されるところである。

以上の構成による解析装置１００の動作を説明する。
図３は、解析装置１００における一連の処理の一例を示すフローチャートである。粒子系取得部１１０は、ＲＭＤ法に倣って繰り込まれた粒子系を取得する（Ｓ２０２）。力演算部１２２は、粒子間の距離から粒子に働く力を演算する（Ｓ２０４）。粒子状態演算部１２４は、演算された力を含む粒子の運動方程式から粒子の速度と位置を演算する（Ｓ２０６）。状態更新部１２６は、粒子データ保持部１５０に保持される粒子の位置、速度を演算された位置、速度で更新する（Ｓ２０８）。終了条件判定部１２８は、終了条件が満たされるか否かを判定する（Ｓ２１０）。終了条件が満たされない場合（Ｓ２１０のＮ）、処理はＳ２０４に戻される。終了条件が満たされる場合（Ｓ２１０のＹ）、表示制御部１３０は演算結果を出力装置１０４に表示させる（Ｓ２１２）。

本実施の形態に係る解析装置１００によると、粒子密度がより低い気相４の粒子に原子径よりも大きい粒子径が設定された粒子系が取得される。大きい粒子径が設定された気相４は連続体として振る舞う。したがって、繰り込み階数を多重に実施しても、液相２と気相４との界面で起こる特徴的な現象を再現できる。
また、本実施の形態に係る解析装置１００によると、気相４の粒子数が増大しないため、計算効率がよい。また、時間刻みを従来法と同程度に維持することが可能であるため、計算効率がよい。また、二相流におけるクヌーセン数を任意に設定できる。

本発明者は、本実施の形態に係る手法の確認計算を行った。図４（ａ）〜（ｔ）は、本実施の形態に係る手法を使用した演算結果を示す図である。図５（ａ）〜（ｔ）は、本実施の形態に係る手法を使用しない場合（従来の手法の場合）の演算結果を示す図である。図４（ａ）〜（ｔ）および図５（ａ）〜（ｔ）は、図１において壁３を取り外した場合の圧力分布の時間変化を示している。

（演算条件）
・気体粒子（本実施の形態に係る手法）
粒子系ｄ：１．２４７２６［ｎｍ］
質量ｍ：１８[ｇ]
ポテンシャルパラメータσ：１．４［ｎｍ］
ポテンシャルパラメータε：２９６［Ｋ］
・気体粒子（従来の手法）
粒子系ｄ：０．４１６５５［ｎｍ］
質量ｍ：１４[ｇ]
ポテンシャルパラメータσ：０．３７１１［ｎｍ］
ポテンシャルパラメータε：７８．６２［Ｋ］
・液体粒子（本実施の形態に係る手法および従来の手法）
粒子系ｄ：０．４［ｎｍ］
質量ｍ：１８[ｇ]
ポテンシャルパラメータσ：０．２６４１［ｎｍ］
ポテンシャルパラメータε：４０４．７１［Ｋ］

本実施の形態に係る手法によると、例えば図４（ｎ）に示されるように、レイリ−テイラー不安定性が現れており、気相の粒子が連続体として振る舞っていることが分かる。これに対して本実施の形態に係る手法を使用しない従来のＲＭＤ法の場合、例えば図５（ｊ）に示されるように、液体粒子が弾道的に散らばっており、気相の粒子が連続体として振る舞っていないことが分かる。

以上、実施の形態に係る解析装置１００の構成と動作について説明した。これらの実施の形態は例示であり、その各構成要素や各処理の組み合わせにいろいろな変形例が可能なこと、またそうした変形例も本発明の範囲にあることは当業者に理解されるところである。

実施の形態では、気液二相流の界面に特徴的な現象をＲＭＤ法により解析する場合について説明したが、これに限られず、様々な混相流の界面に特徴的な現象をＲＭＤ法により解析する場合に適用できる。図６は、変形例に係る解析対象を示す模式図である。図６では、壁としての固相６に気相７が吹き付けている様子を示す。本実施の形態に係る手法は、こうした固気二相流にも適用できる。

１００ 解析装置、 １０２ 入力装置、 １０４ 出力装置、 １１０ 粒子系取得部、 １２０ 数値演算部、 １２２ 力演算部、 １２４ 粒子状態演算部、 １２６ 状態更新部、 １２８ 終了条件判定部、 １３０ 表示制御部、 １５０ 粒子データ保持部。

Claims (4)

1. 繰り込み群分子動力学法を使用して混相流を解析する解析装置であって、
   前記混相流を記述する複数の粒子を含む粒子系を取得する粒子系取得部と、
   粒子間の距離に基づき粒子に働く力を演算する力演算部と、
   離散化された粒子の運動方程式に前記力演算部によって演算された力を適用することによって粒子の位置および速度のうちの少なくともひとつを演算する粒子状態演算部と、を備え、
   前混相流は、少なくとも、第１相と当該第１相より粒子密度が低い第２相とが混在した流れであり、
   前記第１相の粒子には原子径が設定され、前記第２相の粒子には、原子径よりも大きく、クヌーセン数が０．１６以下となり、かつ、粒子同士が重ならないような粒子径が設定されていることを特徴とする解析装置。
2. 前記第２相の粒子には、
   を満たす質量が設定されていることを特徴とする請求項１に記載の解析装置。
3. 前記力演算部は、ポテンシャルエネルギ関数を使用して粒子に働く力を演算する際、前記第２相の粒子に及ぶ引力をゼロとすることを特徴とする請求項１または２に記載の解析装置。
4. 繰り込み群分子動力学法を使用して混相流を解析する機能をコンピュータに実現させる解析コンピュータプログラムであって、
   前記混相流を記述する複数の粒子を含む粒子系を取得する機能と、
   粒子間の距離に基づき粒子に働く力を演算する機能と、
   離散化された粒子の運動方程式に前記演算する機能によって演算された力を適用することによって粒子の位置および速度のうちの少なくともひとつを演算する機能と、を備え、
   前記取得する機能が取得する粒子系に含まれる複数の粒子のうち、密度が高い方の相の粒子には原子径が設定され、密度の低い方の相の粒子には、原子径よりも大きく、クヌーセン数が０．１６以下となり、かつ、粒子同士が重ならないような粒子径が設定されていることを特徴とするコンピュータプログラム。