(57)【要約】【課題】量子流体の粘度に関する情報を得ることができる粘度算出方法を提供する。【解決手段】量子流体の複数の粒子の各々を複数のビーズで表し、複数のビーズの挙動を、経路積分分子動力学法を用いて解析する際に、量子流体に対して、周期境界条件を適用可能な流速分布を生じさせる操作を行う。解析によって得られた量子流体の流速分布と、ポアズイユ流れの流速分布とを対比し、対比結果に基づいて、量子流体の粘度情報を得る。【選択図】図６
2020134144.tif 000002
【特許請求の範囲】
【請求項１】
量子流体の粘度に関する情報を求める方法であって、
前記量子流体の複数の粒子の各々を複数のビーズで表し、前記複数のビーズの挙動を、経路積分分子動力学法を用いて解析する際に、前記量子流体に対して、周期境界条件を適用可能な流速分布を生じさせる操作を行い、
解析によって得られた前記量子流体の流速分布と、ポアズイユ流れの流速分布とを対比し、
対比結果に基づいて、前記量子流体の粘度情報を得る粘度算出方法。
【請求項２】
前記操作は、第１方向に平行な２枚の壁面の間に前記第１方向に平行な流れを生じさせる操作であり、当該操作によって、前記壁面の中央の平面を境界面として一方の空間と他方の空間とに、相互に反平行の向きの流れを生じさせる請求項１に記載の粘度算出方法。
【請求項３】
前記ポアズイユ流れの流速分布として、前記境界面及び前記壁面の位置において流速がゼロになる放物線の分布を採用する請求項２に記載の粘度算出方法。
【請求項４】
量子流体の複数の粒子の各々を複数のビーズで表す条件、解析モデルの境界条件となる平行な２枚の壁面で挟まれた空間を定義する情報を取得する機能と、
前記複数のビーズの挙動を、経路積分分子動力学法を用いて解析する際に、前記量子流体に対して、前記壁面に対して垂直な方向に関して周期境界条件を適用可能な流速分布を生じさせる操作を行う機能と、
解析によって得られた前記量子流体の流速の分布と、ポアズイユ流れの流速分布との対比結果から得られた前記量子流体の粘度情報を出力する機能と
を有する粘度算出装置。
【請求項５】
前記操作は、第１方向に平行な２枚の壁面の間に前記第１方向に平行な流れを生じさせる操作であり、当該操作によって前記壁面の中央の平面を境界面として一方の空間と他方の空間とに、相互に反平行の向きの流れを生じさせる請求項４に記載の粘度算出装置。
【請求項６】
前記ポアズイユ流れの流速分布として、前記境界面及び前記壁面の位置において流速がゼロになる放物線の分布を採用する請求項５に記載の粘度算出装置。
【請求項７】
量子流体の複数の粒子の各々を複数のビーズで表す条件、解析モデルの境界条件となる壁面で挟まれた空間を定義する情報を取得する機能と、
前記複数のビーズの挙動を、経路積分分子動力学法を用いて解析する際に、前記量子流体に対して、前記壁面に対して垂直な方向に関して周期境界条件を適用可能な流速分布を生じさせる操作を行う機能と、
解析によって得られた前記量子流体の流速の分布と、ポアズイユ流れの流速分布との対比結果から得られた前記量子流体の粘度情報を出力する機能と
をコンピュータに実行させるプログラム。
【請求項８】
前記操作は、第１方向に平行な２枚の壁面の間に前記第１方向に平行な流れを生じさせる操作であり、当該操作によって前記壁面の中央の平面を境界面として一方の空間と他方の空間とに、相互に反平行の向きの流れを生じさせる請求項７に記載のプログラム。
【請求項９】
前記ポアズイユ流れの流速分布として、前記境界面及び前記壁面の位置において流速がゼロになる放物線の分布を採用する請求項８に記載のプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、量子流体の粘度算出方法、粘度算出装置、及びこの粘度算出方法をコンピュータに実現させるプログラムに関する。
【背景技術】
【０００２】
分子動力学法を用いたシミュレーションが、現象の微視的な挙動の解析に用いられる。個々の粒子の従う運動方程式を書き下すことができれば、異なる相の連成解析、臨界点近くの様子、乱流等を調べることが可能である。特に、高レイノルズ数の系の解析を行うことが求められている。
【０００３】
レイノルズ数は、系の大きさと流れの速度に比例し、流体の粘度に反比例する。したがって、系を大きくするか、または流体の速度を大きくすることで、高レイノルズ数の系のシミュレーションを行うことが可能になる。ところが、このような方法では、シミュレーション対象の系そのものが大きくなること、及び時間刻み幅を細かくする必要があることから、計算時間が膨大になってしまう。
【０００４】
低粘度流体を表すポテンシャルを見つけることができれば、系の大きさを変えることなく、高レイノルズ数の系を表現することができる。超流動現象は、極低温において流体の粘度がほぼゼロになる現象である。例えば、液体ヘリウム４では、絶対温度２Ｋ程度以下において粘度が急激に低下し、粘度がほぼゼロになる。粘度がほぼゼロになった状態を分子動力学法によって再現することができれば、高レイノルズ数の系のシミュレーションの一助になる。
【０００５】
下記の非特許文献１に、量子系のシミュレーションを行う方法の一つである経路積分分子動力学法（以下、ＰＩＭＤ法という。）について記載されている。また、下記の非特許文献２に、ニュートン流体のせん断粘度をシミュレーションによって求める方法（一様運動量ソース−シンク法（以下、ＵＳＳ法という。））について説明されている。
【先行技術文献】
【非特許文献】
【０００６】
【非特許文献１】S. Miura and S. Okazaki, “Path integral molecular dynamics for bose-einstein and fermi-dirac statistics”, J. Chemical Physics, 112(23) (2000), pp.10116-10124
【非特許文献２】Cao Bingyang and Ruo-Yu Dong, “Nonequilibrium molecular dynamics simulation of shear viscosity by a uniform momentum source-and-sink scheme”, J. Computational Pysics 231(16) (2012), pp.5306-5316
【発明の概要】
【発明が解決しようとする課題】
【０００７】
ＵＳＳ法は、２枚の平行な壁面の間を流れる非圧縮性ニュートン流体の流速分布が放物線になることを利用している。ところが、非ニュートン流体である量子流体の流速分布は放物線にならないため、ＵＳＳ法を用いて粘度を求めることはできない。
【０００８】
本発明の目的は、量子流体の粘度に関する情報を得ることができる粘度算出方法、粘度算出装置、及びプログラムを提供することである。
【課題を解決するための手段】
【０００９】
本発明の一観点によると、
量子流体の粘度に関する情報を求める方法であって、
前記量子流体の複数の粒子の各々を複数のビーズで表し、前記複数のビーズの挙動を、経路積分分子動力学法を用いて解析する際に、前記量子流体に対して、周期境界条件を適用可能な流速分布を生じさせる操作を行い、
解析によって得られた前記量子流体の流速分布と、ポアズイユ流れの流速分布とを対比し、
対比結果に基づいて、前記量子流体の粘度情報を得る粘度算出方法が提供される。
【００１０】
本発明の他の観点によると、
量子流体の複数の粒子の各々を複数のビーズで表す条件、解析モデルの境界条件となる平行な２枚の壁面で挟まれた空間を定義する情報を取得する機能と、
前記複数のビーズの挙動を、経路積分分子動力学法を用いて解析する際に、前記量子流体に対して、前記壁面に対して垂直な方向に関して周期境界条件を適用可能な流速分布を生じさせる操作を行う機能と、
解析によって得られた前記量子流体の流速の分布と、ポアズイユ流れの流速分布との対比結果から得られた前記量子流体の粘度情報を出力する機能と
を有する粘度算出装置が提供される。
【００１１】
本発明のさらに他の観点によると、
量子流体の複数の粒子の各々を複数のビーズで表す条件、解析モデルの境界条件となる壁面で挟まれた空間を定義する情報を取得する機能と、
前記複数のビーズの挙動を、経路積分分子動力学法を用いて解析する際に、前記量子流体に対して、前記壁面に対して垂直な方向に関して周期境界条件を適用可能な流速分布を生じさせる操作を行う機能と、
解析によって得られた前記量子流体の流速の分布と、ポアズイユ流れの流速分布との対比結果から得られた前記量子流体の粘度情報を出力する機能と
をコンピュータに実行させるプログラムが提供される。
【発明の効果】
【００１２】
上述の方法により、量子流体の粘度情報を得ることができる。
【図面の簡単な説明】
【００１３】
【図１】図１は、ＨＦＤ−Ｂ２ポテンシャルの形状を示すグラフである。
【図２】図２Ａ及び図２Ｂは、それぞれ温度が４Ｋ及び０．５Ｋのときのビーズの位置のＰＩＭＤ法による解析結果を示す斜視図である。
【図３】図３Ａは、ＵＳＳ法で用いられる解析モデルを示す図であり、図３Ｂは、ＵＳＳ法を適用して得られる流速分布を説明するための図である。
【図４】図４Ａは、常流体の流速分布を示すグラフであり、図４Ｂは、超流体の流速分布を示すグラフであり、図４Ｃは、図４Ａと図４Ｂとのグラフを合わせた超流動状態の流体の流速分布を示すグラフである。
【図５】図５は、ｘ＝０及びｘ＝ａの位置に、ｙｚ面に平行に配置された２枚の壁面の間の量子流体の流速分布と、ポアズイユ流れの流速分布とを重ねて示すグラフである。
【図６】図６は、実施例による粘度算出方法のフローチャートである。
【図７】図７Ａ及び図７Ｂは、ＰＩＭＤ法及び古典的分子動力学法にＵＳＳ法を適用して、量子流体及びニュートン流体の流速分布を求めた結果を示すグラフである。
【図８】図８は、量子流体の温度ごとに実施例による方法で求められた粘度比ηを示すグラフである。
【図９】図９は、実施例による粘度算出装置のブロック図である。
【発明を実施するための形態】
【００１４】
実施例について説明する前に、実施例で用いるＰＩＭＤ法及びＵＳＳ法について説明する。ＰＩＭＤ法では、量子的な１個の粒子をＰ個のビーズ（レプリカ）で表現し、これらのビーズの挙動を、古典的な分子動力学法によって求める。
【００１５】
［ＰＩＭＤ法］
まず、量子的な系の１個の粒子に対応するＰ個のビーズの古典系のハミルトニアンについて説明する。
【００１６】
量子的な系の分配関数ＺＱＭは、以下の式で表される。
【数１】
2020134144.tif 000003
ここで、Ｈハットはハミルトニアンであり、Ｔは温度であり、ｋＢはボルツマン定数である。記号Ｔｒは、行列の対角和を意味する。
【００１７】
量子的な１個の粒子は、古典系ではＰ個のビーズによって表現される。ｓ番目のビーズの座標をｒ（ｓ）（ｓ＝１、２、・・・Ｐ）と表記する。便宜上、ｒ（Ｐ＋１）＝ｒ（１）とする。ポテンシャルＶ（ｒ）中の１粒子系では、古典的分配関数ＺＰは以下の式で表される。
【数２】
2020134144.tif 000004
ここで、Ｐは整数、ｍは量子的な粒子の質量を表す。１つの粒子に対応するＰ個のビーズは、有効ポテンシャルＶｅｆｆを表す式の右辺の第二項に示すように、バネによって環状につながっている。
【００１８】
式（１）に示した量子系の分配関数ＺＱＭは、式（２）に示した古典的分配関数ＺＰにおいてパラメータＰが無限大の極限を取ったものと一致する。すなわち、以下の関係が成立する。
【数３】
2020134144.tif 000005
ＰＩＭＤ法では、古典的分配関数ＺＰのパラメータＰを有限の値に設定して古典的運動方程式を解く。
【００１９】
上述の式（２）の右辺には、系のポテンシャルのみが明示的に表れている。ｓ番目のビーズの仮想的運動量ｐ（ｓ）と各ビーズの質量Ｍを導入して、式（２）は以下のように書き換えることができる。
【数４】
2020134144.tif 000006
ここで、Ｈｅｆｆは、以下の式（５）で表される。
【数５】
2020134144.tif 000007
【００２０】
式（５）で表されたＨｅｆｆは、元の量子系に対応した古典系のハミルトニアンとみなすことができる。式（５）の右辺の第１項は、ビーズの運動エネルギを表し、第２項は元の粒子間に作用するポテンシャルエネルギをビーズ間に作用するように置き換えたポテンシャルエネルギを表し、第３項は、１つの粒子に対応するＰ個のビーズの間に作用するバネ力（引力）によるポテンシャルエネルギを表している。このように、１つの粒子に対応するＰ個のビーズの間に、元の粒子の間に作用するポテンシャルＶ（ｒ）に基づいて得られるポテンシャル（式（５）の右辺第二項）と、ビーズの間に引力を生じさせるポテンシャル（式（５）の右辺第三項）とを重畳した相互作用ポテンシャルが作用する。
【００２１】
以下に説明する実施例では、Ｍ＝ｍ／Ｐとする。すなわち、Ｐ個のビーズの質量の合計が、量子的な１個の粒子の質量に等しい。
【００２２】
量子的な粒子が複数の場合には、粒子の統計性を考慮しなければならない。同種の粒子がＮ個である場合、その状態空間を張る完全系は、以下のように表される。
【数６】
2020134144.tif 000008
ここで、ｒは粒子の位置であり、Ｓ-ＮはＮ次対称群であり、σは置換である。量子的な粒子がボソンのときε＝１である。量子的な粒子がフェルミオンのとき、σが偶ならε＝１、σが奇ならε＝−１である。
【００２３】
このとき、式（５）の最近接バネの部分が変更を受ける。ｉ番目の量子的な粒子のｓ番目のビーズをｒｉ（ｓ）（ｉ＝１，２，・・・，Ｎ、ｓ＝１，２，・・・，Ｐ）と表記し、以下の行列要素Ａｉｊ（ｓ）を定義する。
【数７】
2020134144.tif 000009
【００２４】
量子的な系の分配関数ＺＱＭは、以下のように計算することができる。
【数８】
2020134144.tif 000010
【００２５】
ここで、Ｄは空間の次元数である。ＺＰは、古典的な系の分配関数であり、以下の式で表される。
【数９】
2020134144.tif 000011
式（９）のうち、
【数１０】
2020134144.tif 000012
の項は、量子的な粒子がフェルミオンのときは行列式ｄｅｔＡｉｊ（ｓ）であり、ボソンのときはｐｅｒｍＡｉｊ（ｓ）である。これを用いると、古典的な分配関数ＺＰは以下のように書くことができる。
【数１１】
2020134144.tif 000013
【００２６】
量子的な粒子がフェルミオンのとき、
【数１２】
2020134144.tif 000014
と表され、量子的な粒子がボソンのとき、
【数１３】
2020134144.tif 000015
と表される。
【００２７】
式（１２）、式（１３）が、多粒子系でのＰＩＭＤ法で用いる有効ポテンシャルである。一粒子系においては、Ｐ個のビーズが環状に連なっているが、多粒子系においては、そのような単純な形にはならない。有効ポテンシャルはｄｅｔＡｉｊ（ｓ）またはｐｅｒｍＡｉｊ（ｓ）を含んでおり、この計算にはΟ（Ｎ！）程度の時間を要する。ここで、Οはランダウの漸近記号である。
【００２８】
次に、図１〜図２Ｂを参照して、ＰＩＭＤ法を用いてヘリウム４の粒子の挙動をシミュレーションした結果について説明する。
【００２９】
ヘリウム４の原子間ポテンシャルとして、ＨＦＤ−Ｂ２ポテンシャルを用いる。これは、古典的液体へリウムを記述する有効ポテンシャルである。ＨＦＤ−Ｂ２ポテンシャルＶＨＦＤＢ２は、以下の式で定義される。なお、以下の式に現れるＤは、式（１１）の空間の次元数Ｄとは異なる。
【数１４】
2020134144.tif 000016
【００３０】
図１は、ＨＦＤ−Ｂ２ポテンシャルの形状を示すグラフである。横軸は距離ｒを単位「ｎｍ」で表し、縦軸はＨＦＤ−Ｂ２ポテンシャルＶＨＦＤＢ２をボルツマン定数ｋＢで除した値を単位（Ｋ）で表す。ＨＦＤ−Ｂ２ポテンシャルは、分子間ポテンシャルとしてしばしば用いられるレナードジョーンズポテンシャルと似た形状を持つ。距離ｒ＝０．２５ｎｍの近傍に、強い斥力による壁が存在し、距離ｒ＝０．３ｎｍの近傍にポテンシャルの底が存在する。ポテンシャルが底を示す距離より遠い距離では引力が働き、距離ｒ→∞でポテンシャルがゼロになる。
【００３１】
式（１３）の右辺第１項にＨＦＤ−Ｂ２ポテンシャルを適用したボソンの有効ポテンシャルを用いて、ヘリウム４の原子の挙動をＰＩＭＤ法により解析した。原子数Ｎ＝７とし、ビーズ数Ｐ＝１２８とした。
【００３２】
図２Ａ及び図２Ｂは、それぞれ温度が４Ｋ及び０．５Ｋのときのビーズの位置の解析結果を示す斜視図である。温度が４Ｋのときは複数のビーズがある程度の塊となっており、個々の原子を認識することができる。ところが、温度が０．５Ｋのときには、複数のビーズはばらばらに分布しており、個々の原子を認識することはできない。
【００３３】
［ＵＳＳ法］
次に、図３Ａ及び図３Ｂを参照して、通常の古典的分子動力学法で人工的にポアズイユ流れを作り出し、その流速分布から粘度を計算するＵＳＳ法について説明する。
【００３４】
図３Ａは、ＵＳＳ法で用いられる解析モデルを示す図である。ｙｚ平面と平行に置かれた２枚の壁面１０の間を、ｚ軸方向に一様に流れる非圧縮性のニュートン流体の定常状態を考える。平行に配置された２枚の壁面１０の間に流体の複数の粒子１１が配置される。２枚の壁面１０の一方のｘ座標をｘ＝０、他方のｘ座標をｘ＝ａとする。すなわち、２枚の壁面１０の間隔はａである。流束ｖと圧力ｐとは、以下の式を満たす。
【数１５】
2020134144.tif 000017
【００３５】
このとき、ナビエストークス方程式は、粘度をμとして以下の式で表される。
【数１６】
2020134144.tif 000018
【００３６】
ｘ＝０及びｘ＝ａの位置においてｖｚ＝０となる境界条件を適用すると、式（１５）、式（１６）の解は、
【数１７】
2020134144.tif 000019
となる。これは、ポアズイユ流れと呼ばれる。式（１７）は、ｚ軸方向に圧力勾配が存在する場合に、放物線型の流速分布が生じることを意味する。圧力勾配と流速分布がわかれば、粘度μを求めることができる。
【００３７】
分子動力学法では、粒子の速度等に周期境界条件を適用することが便利である。ＵＳＳ法では、上述の境界条件を満たし、かつ周期境界条件を満たすような流速分布を生じさせる。
【００３８】
図３Ｂは、ＵＳＳ法を適用して得られる流速分布を説明するための図である。２枚の壁面１０で挟まれた空間を、ｙｚ面に平行な境界面１２で二等分する。二等分された２つの領域中の粒子１１に対して、一定時間ごとに以下の第１〜第３操作を実行する。
【００３９】
第１操作では、一方の領域（右半分の領域）から速度のｚ成分が最大の粒子１１を選択する。第２操作では、他方の領域（左半分の領域）から速度のｚ成分ｖｚが最小（ｖｚが負で絶対値が最大）の粒子１１を選択する。第３操作では、選択された２つの粒子１１の速度のｚ成分ｖｚの値を交換する。このとき、速度のｘ成分ｖｘ、ｙ成分ｖｙ、及び粒子の座標の値は交換しない。
【００４０】
第１操作〜第３操作は、一定時間ごとに、右半分の領域の粒子１１にｚ軸の負の方向に撃力を作用させ、左半分の領域の粒子１１には、ｚ軸の正の方向に撃力を作用させていると言い換えることができる。すなわち、右半分の領域と左半分の領域とでは、ｚ軸方向に関して相互に反対向きの圧力勾配が生じる。従って、解析対象の流体がニュートン流体であれば、定常状態において右半分と左半分とで、ｚ軸方向に関して相互に反対向きの放物線型の流速分布１３が生じる。このため、ｘ軸方向に関して、壁面１０の位置において周期境界条件を適用することができる。なお、実際の計算では、ｙ軸方向及びｚ軸方向についても周期境界条件を適用する。例えば、解析領域として、ｘ軸方向、ｙ軸方向、ｚ軸方向に垂直な面を持つ立方体を採用するとよい。計算で得られた流速分布１３に基づいて、式（１７）に従って粘度μを求めることができる。
【００４１】
［実施例］
次に、実施例によるシミュレーション方法について説明する。
ＰＩＭＤ法で作り出した超流動状態の流体にＵＳＳ法を適用して量子流体の流れを生じさせ、量子流体の粘度を求めたい。ところが、超流動状態の流体は非ニュートン流体であり、流速分布は放物線型にならないため、そのままＵＳＳ法を適用して粘度を求めることはできない。本実施例では、超流動状態をよく記述する二流体モデルを利用する。
【００４２】
二流体モデルでは、超流動状態の液体ヘリウム４を、ゼロでない通常の粘度を持つ常流体（常流動成分）と、粘度がゼロである超流体（超流動成分）との混合物として記述する。常流体の速度場をｖｎ、粒子数比をｎｎと表し、超流体の速度場をｖｓ、粒子数比をｎｓと表す。このとき、流体全体の速度場ｖは以下の式で表される。
【数１８】
2020134144.tif 000020
【００４３】
液体ヘリウム４は、臨界温度２．１７Ｋより高温では常流動であるため、臨界温度２．１７Ｋより高温のとき粒子数比は以下の値をとる。
【数１９】
2020134144.tif 000021
【００４４】
逆に、液体ヘリウム４の温度が０Ｋに近づくと、粒子数比は以下の値に漸近する。
【数２０】
2020134144.tif 000022
【００４５】
この二流体モデルにＵＳＳ法を適用した場合、流速分布がどのようになるか検討する。常流体は通常の粘性流体であるから、常流体の流速ｖｎの分布は放物線を描くと期待される。一方、超流体は粘性を持たないので、流速分布は平坦である。常流体と超流体とを合わせた流体の流速分布は、放物線型の流速分布と平坦な流速分布とを合わせた形状になると期待される。
【００４６】
図４Ａは、常流体の流速分布を示すグラフであり、図４Ｂは、超流体の流速分布を示すグラフである。図４Ｃは、図４Ａと図４Ｂとのグラフを合わせた超流動状態の流体の流速分布を示すグラフである。二流体モデルの粘性は常流体の占める割合そのものであるから、この系の粘性は、流速分布が急激に変化する位置から読み取ることができる。
【００４７】
次に、図５を参照して、量子流体の粘度を計算する方法について説明する。
図５は、ｘ＝０及びｘ＝ａの位置に、ｙｚ面に平行に配置された２枚の壁面１０（図３Ｂ）の間の量子流体の流速分布２１と、ポアズイユ流れの流速分布２２とを重ねて示すグラフである。ＰＩＭＤ法による解析にＵＳＳ法を適用することにより、量子流体の流速分布２１が得られる。ポアズイユ流れの流速分布２１は、古典的分子動力学法にＵＳＳ法を適用して得られたものでもよいし、理論式から得られたものでもよい。ポアズイユ流れの流速分布２１は放物線型になる。
【００４８】
ＰＩＭＤ法による解析にＵＳＳ法を適用すると、０＜ｘ＜（ａ／２）の範囲で流速ｖが正になり、（ａ／２）＜ｘ＜ａの範囲で流速ｖが負になる流速分布２１が求まる。流速分布２１は、ｘ＝（ａ／２）の位置において急激に変化する。放物線型の流速分布２２は、その最大値が流速分布２１の最大値と等しくなるように決定されている。
【００４９】
あるｘ座標ｘ０における流速分布２１と最大値との差をｄ１、流速分布２２と最大値との差をｄ２で表す。量子流体の粘度情報として、量子流体の全体に対する常流体の成分の比（以下、粘度比ηという。）を採用する。粘度比ηは、以下の式で表すことができる。
【数２１】
2020134144.tif 000023
【００５０】
ｘ軸上の複数の位置で粘度比ηを計算し、その平均値を量子流体の粘度比ηとするとよい。なお、ｘ＝０、ａ／２、ａの近傍、及び流速が最大値を取る位置、すなわちｘ＝ａ／４、３ａ／４の近傍では揺らぎによって粘度比ηの値が不安定になるため、これらの座標の近傍の値は、粘度比ηの平均値の計算対象から除外するとよい。
【００５１】
次に、図６〜図７Ｂを参照して、実施例による粘度算出方法について説明する。
図６は、実施例による粘度算出方法のフローチャートである。まず、粘度を算出する量子流体の粒子を定義する情報、境界条件、初期条件等を決定する（ステップＳ１）。量子流体の粒子を定義する情報には、量子的な粒子の質量ｍ及び個数Ｎ、量子的な１個の粒子に対応するビーズの個数Ｐ、粒子によるポテンシャル（式（１２）、式（１３））を定義する情報等が含まれる。
【００５２】
解析領域内に複数の粒子に対応する複数のビーズを配置し、ＰＩＭＤ法により、ビーズの位置及び速度を、１タイムステップ分だけ時間発展させる（ステップＳ２）。なお、一定時間ごとに温度制御を行う。温度制御には、例えば、能勢フーバーチェイン法を用いる。現時点のタイムステップで、ＵＳＳ法による操作、すなわち周期境界条件を適用可能な流速分布を生じさせる操作を実行するか否かを判定する（ステップＳ３）。例えば、一定のタイムステップ数ごとにＵＳＳ法による操作を実行する。ＵＳＳ法による操作を実行する場合には、図３Ｂを参照して説明した反平行の流れを生じさせる操作、すなわち粒子の速度のｚ成分の値を交換する操作を実行する（ステップＳ４）。この操作を行った後、次のタイムステップの時間発展を行う（ステップＳ２）。ＵＳＳ法による操作を実行しない場合には、流体の流れが定常状態に到達したか否かを判定する（ステップＳ５）。定常状態に到達していない場合には、次のタイムステップの時間発展を行う（ステップＳ２）。
【００５３】
流体の流れが定常状態に到達していると判定されたら、ステップＳ２で行っていた温度制御を停止する（ステップＳ６）。その後、ＰＩＭＤ法によるビーズの位置及び速度の時間発展（ステップＳ７）、ＵＳＳ法による操作を実行するか否かの判定（ステップＳ８）を行う。ステップＳ７では、温度制御は行わない。ＵＳＳ法による操作を実行しない場合には、流速分布のサンプリングを行う（ステップＳ１０）。ＵＳＳ法による操作を実行する場合には、ステップＳ４と同様に反平行の流れを生じさせる操作を実行（ステップＳ９）した後、流速分布のサンプリングを行う（ステップＳ１０）。
【００５４】
流速分布のサンプリングを行った後、サンプル数が十分であるか否かを判定する（ステップＳ１１）。サンプル数が不十分である場合には、ＰＩＭＤ法による次のタイムステップの時間発展を行う（ステップＳ７）。サンプル数が十分である場合には、サンプリングされた複数の状態基づいて、ビーズの速度のｚ成分の値から、ｘ軸方向に関する流速分布を計算する（ステップＳ１２）。次に、サンプリングされた複数の状態の平均の流速分布に基づいて、図５を参照して説明した方法により、量子流体の粘度比ηを計算する（ステップＳ１３）。
【００５５】
次に、実施例による方法で、実際に液体ヘリウム４の粘度を算出した結果について説明する。
【００５６】
図７Ａ及び図７Ｂは、ＰＩＭＤ法及び古典的分子動力学法にＵＳＳ法を適用して、量子流体及びニュートン流体の流速分布を求めた結果を示すグラフである。図７Ａは、量子的な１個の粒子に対応させるビーズの個数Ｐを１６とした場合、及び図７ＢはＰを１２８とした場合の結果を示す。グラフ中の丸記号、四角記号、三角記号、菱形記号、及び星形記号は、それぞれ温度Ｔを１．０Ｋ、１．５Ｋ，２．０Ｋ、２．５Ｋ、１５．０Ｋとした場合の速度を示す。温度Ｔが１５．０Ｋの場合は超流動状態が発現しないため、古典的に分子動力学法を用いて速度を求めた。
【００５７】
ＰＩＭＤ法及び古典的分子動力学法に共通のシミュレーション条件は、以下の通りである。
・ヘリウム４原子の質量：６．６４×１０−２７ｋｇ
・時間刻み：３．０×１０−１６ｓ
・境界条件：周期境界条件
【００５８】
ＰＩＭＤ法のシミュレーション条件は、以下の通りである。
・ヘリウム４原子数：Ｎ＝８
・解析領域：一辺の長さが７．１５９×１０−１０ｍの立方体
・ビーズ数：Ｐ＝１６及び１２８
・ＵＳＳ法による粒子の速度の値の交換頻度
Ｐ＝１２８、Ｔ＝２．５Ｋのとき２９０タイムステップごと
Ｐ＝１２８、Ｔ＝２．０Ｋのとき２６０タイムステップごと
Ｐ＝１２８、Ｔ＝１．５Ｋのとき２００タイムステップごと
Ｐ＝１２８、Ｔ＝１．０Ｋのとき８５タイムステップごと
Ｐ＝１６、Ｔ＝２．５Ｋのとき６２０タイムステップごと
Ｐ＝１６、Ｔ＝２．０Ｋのとき５２０タイムステップごと
Ｐ＝１６、Ｔ＝１．５Ｋのとき５３０タイムステップごと
Ｐ＝１６、Ｔ＝１．０Ｋのとき３８０タイムステップごと
計算条件ごとに速度の値の交換頻度を変えているのは、計算の安定性のためである。
【００５９】
古典的分子動力学法のシミュレーション条件は、以下の通りである。
・原子数：Ｎ＝４００
・解析領域：一辺の長さが２．６３７×１０−９ｍの立方体
・ＵＳＳ法による粒子の速度の値の交換頻度：１０００タイムステップごと
【００６０】
解析領域をｘ方向に６０等分して６０個の区画に区分し、これらの区画に０から５９まで通し番号を付す。区画ごとにビーズのｚ方向の速度の平均値を求め、この平均値を当該区画における量子流体の流速とした。図７Ａ及び図７Ｂの横軸に付した数値は、区画に付された通し番号を表している。図７Ａ、図７Ｂの横軸の−０．５、５９．５の位置が、それぞれ図３Ｂのｘ＝０、ｘ＝ａの位置の壁面１０に相当し、横軸の２９．５の位置が、図３Ｂのｘ＝ａ／２の位置の境界面１２に相当する。
【００６１】
１４番目の区画の流速と１５番目の区画の流速との平均値で、各区画の流速を規格化した。図７Ａ及び図７Ｂの縦軸は、規格化したｚ方向の流速を表す。
【００６２】
図７Ａ及び図７Ｂに示すように、ＰＩＭＤ法にＵＳＳ法を適用して求めた流速分布は、古典的分子動力学法にＵＳＳ法を適用して求めた放物線型の流速分布からずれていることがわかる。
【００６３】
図７Ａ及び図７Ｂに示した解析結果に基づいて、図５を参照して説明した方法を用いて量子流体の粘度比ηを求めた。
【００６４】
図８は、量子流体の温度と、実施例による方法で求められた粘度比ηとの関係を示すグラフである。横軸は温度を単位「Ｋ」で表し、縦軸は粘度比ηを表す。グラフ中の四角記号及び丸記号は、それぞれ量子流体の１個の粒子に対応するビーズの個数Ｐが１６個及び１２８個の場合の計算結果を示す。比較のために、実際の実験によって求められている液体ヘリウム４の粘度比を破線で示す。
【００６５】
実施例による方法で求めた粘度比ηは、低温になるほど低下する傾向を示しており、破線で示した実験結果がある程度反映されていることがわかる。また、ビーズ数Ｐを大きくした方が、温度１Ｋ程度の極低温における粘度比が低くなっている。このように、ビーズ数Ｐを大きくした方が、より高精度に粘度比ηを求めることができる。
【００６６】
次に、図９を参照して実施例による粘度算出装置について説明する。
図９は、実施例による粘度算出装置のブロック図である。粘度算出装置は一般的なコンピュータで構成され、中央処理ユニット（ＣＰＵ、処理装置）５０、メモリ５１、入出力装置５２、及び補助記憶装置５３を含む。これらは、バス５４で相互に接続されている。
【００６７】
入出力装置５２は、キーボード、マウス等のポインティングデバイス、ディスプレイ、リムーバブルメディアのリーダライタ、通信装置等を含む。ディスプレイは、ユーザの操作に必要な各種ウィンドウやデータ等を表示する。通信装置は、外部機器とデータ通信を行う。ユーザは、入出力装置を操作することにより、コンピュータに対する指示、各処理に必要なデータの入力等を行う。また、処理結果が、入出力装置に出力される。
【００６８】
補助記憶装置５３は、ハードディスク等で構成され、コンピュータのＯＳ、粘度算出プログラム、粘度の算出に必要なデータや処理結果等を蓄積する。
【００６９】
メモリ５１は、ＲＯＭ（Ｒｅａｄ  Ｏｎｌｙ  Ｍｅｍｏｒｙ）やＲＡＭ（Ｒａｎｄｏｍ  Ａｃｃｅｓｓ  Ｍｅｍｏｒｙ）等で構成される。メモリ５１は、ＣＰＵ５０により補助記憶装置５３から読み出された粘度算出プログラム等を格納する。
【００７０】
ＣＰＵ５０は、コンピュータのＯＳ、メモリ５１に格納されている粘度算出プログラム等に基づいて、各種演算や入出力装置５２及び補助記憶装置５３の制御を行う。
【００７１】
次に、粘度算出装置の動作について説明する。ステップＳ１（図６）で決定された粒子の定義、境界条件、初期条件等が入出力装置５２に入力され、ＣＰＵ５０がこれらの情報を取得する。その後、ＣＰＵ５０は、ステップＳ２からステップＳ７までの処理を実行する。その後、ステップＳ７で求められた粘度の値をディスプレイに表示させる。
【００７２】
次に、上記実施例の優れた効果について説明する。
従来は、古典的分子動力学法による解析にＵＳＳ法を適用して流体の粘度情報を求めていたが、超流動状態を示す量子流体の粘度情報を求めることはできなかった。上記実施例の方法を用いることにより、超流動状態を示す量子流体の粘度情報を求めることが可能になる。
【００７３】
上述の実施例は例示であり、本発明は上述の実施例に制限されるものではない。例えば、種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。
【符号の説明】
【００７４】
１０  壁面
１１  粒子
１２  境界面
１３、２１、２２  流速分布
５０  ＣＰＵ
５１  メモリ
５２  入出力装置
５３  補助記憶装置
５４  バス
【図１】
2020134144.tif 000024
【図２】
2020134144.tif 000025
【図３】
2020134144.tif 000026
【図４】
2020134144.tif 000027
【図５】
2020134144.tif 000028
【図６】
2020134144.tif 000029
【図７】
2020134144.tif 000030
【図８】
2020134144.tif 000031
【図９】
2020134144.tif 000032