(57)【要約】【課題】自由表面を有する流体等のシミュレーションを行うことが可能なシミュレーション方法を提供する。【解決手段】粒子系に含まれる粒子の挙動を、各粒子について運動方程式を数値的に解くことにより粒子系の状態を変化させるタイムステップを複数回繰り返す。まず、粒子密度の設定値、及び設定温度を決定する。タイムステップにおいて、粒子系の状態に基づいて、粒子の各々の位置を含む局所領域の粒子密度を算出し、粒子密度の算出値と粒子密度の設定値とに基づいて、設定温度を更新し、更新された設定温度に基づいて温度制御を行う。【選択図】図１
2019125296.tif 000002
【特許請求の範囲】
【請求項１】
粒子系に含まれる複数の粒子の各々について運動方程式を数値的に解くことにより前記粒子系の状態を変化させるタイムステップを複数回繰り返すシミュレーション方法であって、
粒子密度の設定値、及び設定温度を決定し、
前記タイムステップにおいて、
前記粒子系の状態に基づいて、前記粒子の各々の位置を含む局所領域の粒子密度を算出し、
粒子密度の算出値と粒子密度の前記設定値とに基づいて、前記設定温度を更新し、
更新された前記設定温度に基づいて温度制御を行うシミュレーション方法。
【請求項２】
前記粒子に重力が作用しており、前記粒子系に特性速度が定義され、
前記粒子について運動方程式を解く際に、前記粒子に作用させる重力を実際の重力より大きくし、かつ運動方程式を解く際に用いる特性速度を実際の特性速度より速くする請求項１に記載のシミュレーション方法。
【請求項３】
前記粒子に作用させる重力を実際の重力より大きくし、運動方程式を解く際に用いる特性速度を実際の特性速度より速くする際に、前記粒子系のフルード数が変化しないように、重力及び特性速度を変化させる請求項２に記載のシミュレーション方法。
【請求項４】
前記設定温度を更新する処理において、粒子密度の前記算出値が前記設定値より高いとき、前記設定温度を現在の値より高くし、粒子密度の前記算出値が前記設定値より低いとき、前記設定温度を現在の値より低くする請求項１乃至３のいずれか１項に記載のシミュレーション方法。
【請求項５】
前記設定温度の加速度項と、粒子密度の前記算出値と前記設定値との差に比例する外力と、前記粒子の質量に比例する粘性摩擦力を含む運動方程式を解くことにより、前記設定温度を更新する請求項１乃至４のいずれか１項に記載のシミュレーション方法。
【請求項６】
解析領域に複数の粒子を配置して、複数の前記粒子からなる粒子系の挙動を、分子動力学法によりシミュレーションする処理手段を有するシミュレーション装置であって、
前記処理手段は、
入力手段から、粒子密度の設定値を特定するパラメータ及び設定温度を特定するパラメータを含むシミュレーションの初期設定条件を取得する機能、
入力された前記初期設定条件に基づいて、前記粒子系の挙動をシミュレーションする機能、
前記粒子の各々の位置を含む局所領域の粒子密度を算出する機能、
粒子密度の算出値と粒子密度の前記設定値とに基づいて、前記設定温度を更新する機能、及び
更新された前記設定温度に基づいて温度制御を行う機能
を有するシミュレーション装置。
【請求項７】
前記処理手段は、粒子密度の前記算出値が前記設定値より高いとき、前記設定温度を現在の値より高くし、粒子密度の前記算出値が前記設定値より低いとき、前記設定温度を現在の値より低くすることにより、前記設定温度を更新する機能を実現する請求項６に記載のシミュレーション装置。
【請求項８】
前記処理手段は、前記設定温度の加速度項と、粒子密度の前記算出値と前記設定値との差に比例する外力と、前記粒子の質量に比例する粘性摩擦力を含む運動方程式を解くことにより、前記設定温度を更新する機能を実現する請求項６または７に記載のシミュレーション装置。
【請求項９】
解析領域に複数の粒子を配置して、複数の前記粒子からなる粒子系の挙動を、分子動力学法によりシミュレーションする機能をコンピュータに実現させるプログラムであって、
入力手段から、粒子密度の設定値及び設定温度を含むシミュレーションの初期設定条件を取得する機能、
入力された前記初期設定条件に基づいて、前記粒子系の挙動をシミュレーションする機能、
前記粒子の各々の位置を含む局所領域の粒子密度を算出する機能、
粒子密度の算出値と粒子密度の前記設定値とに基づいて、前記設定温度を更新する機能、及び
更新された前記設定温度に基づいて温度制御を行う機能
を有するプログラム。
【請求項１０】
粒子密度の前記算出値が前記設定値より高いとき、前記設定温度を現在の値より高くし、粒子密度の前記算出値が前記設定値より低いとき、前記設定温度を現在の値より低くすることにより、前記設定温度を更新する機能を実現する請求項９に記載のプログラム。
【請求項１１】
前記設定温度の加速度項と、粒子密度の前記算出値と前記設定値との差に比例する外力と、前記粒子の質量に比例する粘性摩擦力を含む運動方程式を解くことにより、前記設定温度を更新する機能を実現する請求項９または１０に記載のプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、シミュレーション方法、シミュレーション装置、及びプログラムに関する。
【背景技術】
【０００２】
分子動力学法に温度制御を取り入れてシミュレーションを行う方法が公知である（例えば、特許文献１）。例えば、シミュレーション対象の粒子系の温度が設定温度に近づくように、分子動力学法によって算出された粒子の速度を、所望の設定温度、及び粒子系の温度に基づいて更新する。
【０００３】
また、従来、分子動力学法ではミクロな領域の解析が主であったため、解析領域の境界条件として周期境界条件が適用されることが多い。周期境界条件を適用すると、解析領域内の粒子密度は不変であるため、粒子密度を制御する手法は採用されていなかった。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特開平１０−３３４０７６号公報
【発明の概要】
【発明が解決しようとする課題】
【０００５】
くりこみ群分子動力学法ではマクロな解析を実施可能であるため、自由表面を有する流体の挙動を解析することも可能になる。自由表面を有する流体が定常状態に達するまでの時間を短くするために、重力加速度及び特性速度を大きくしてシミュレーションを行う場合がある。一例として、無次元数であるフルード数が変化しないように、重力加速度と特性速度とを大きくする。
【０００６】
本願の発明者らが行った重力加速度を大きくした従来のシミュレーションによると、自由表面を有する流体の挙動が、実際の流体の挙動からかけ離れたものになってしまうことが判明した。本発明の目的は、自由表面を有する流体等のシミュレーションを行うことが可能なシミュレーション方法、シミュレーション装置、及びプログラムを提供することである。
【課題を解決するための手段】
【０００７】
本発明の一観点によると、
粒子系に含まれる粒子の挙動を、各粒子について運動方程式を数値的に解くことにより前記粒子系の状態を変化させるタイムステップを複数回繰り返すシミュレーション方法であって、
粒子密度の設定値、及び設定温度を決定し、
前記タイムステップにおいて、
前記粒子系の状態に基づいて、前記粒子の各々の位置を含む局所領域の粒子密度を算出し、
粒子密度の算出値と粒子密度の前記設定値とに基づいて、前記設定温度を更新し、
更新された前記設定温度に基づいて温度制御を行うシミュレーション方法が提供される。
【０００８】
本発明の他の方法によると、
解析領域に複数の粒子を配置して、複数の前記粒子からなる粒子系の挙動を、分子動力学法によりシミュレーションする処理手段を有するシミュレーション装置であって、
前記処理手段は、
入力手段から、粒子密度の設定値を特定するパラメータ及び設定温度を特定するパラメータを含むシミュレーションの初期設定条件を取得する機能、
入力された前記初期設定条件に基づいて、前記粒子系の挙動をシミュレーションする機能、
前記粒子の各々の位置を含む局所領域の粒子密度を算出する機能、
粒子密度の算出値と粒子密度の前記設定値とに基づいて、前記設定温度を更新する機能、及び
更新された前記設定温度に基づいて温度制御を行う機能
を有するシミュレーション装置が提供される。
【０００９】
本発明のさらに他の観点によると、
解析領域に複数の粒子を配置して、複数の前記粒子からなる粒子系の挙動を、分子動力学法によりシミュレーションする機能をコンピュータに実現させるプログラムであって、
入力手段から、粒子密度の設定値及び設定温度を含むシミュレーションの初期設定条件を取得する機能、
入力された前記初期設定条件に基づいて、前記粒子系の挙動をシミュレーションする機能、
前記粒子の各々の位置を含む局所領域の粒子密度を算出する機能、
粒子密度の算出値と粒子密度の前記設定値とに基づいて、前記設定温度を更新する機能、及び
更新された前記設定温度に基づいて温度制御を行う機能
を有するプログラムが提供される。
【発明の効果】
【００１０】
自由表面を有する流体の挙動を、十分な精度でシミュレーションすることが可能になる。
【図面の簡単な説明】
【００１１】
【図１】図１は、実施例によるシミュレーション方法のフローチャートである。
【図２】図２は、実施例による方法でシミュレーションを行ったシミュレーション結果を示す図である。
【図３】図３は、レイノルズ数Ｒｅと動力数Ｎｐとの関係を示すグラフである。
【図４】図４は、実施例によるシミュレーション装置のブロック図である。
【図５】図５Ａは、シミュレーション対象の撹拌槽モデルの縦断面図であり、図５Ｂは図５Ａの一点鎖線５Ｂ−５Ｂにおける横断面図である。
【図６】図６は、従来の分子動力学法を用いて流体の挙動をシミュレーションした結果を示す図である。
【発明を実施するための形態】
【００１２】
実施例について説明する前に、自由表面を有する流体の挙動を従来のシミュレーション方法でシミュレーションした結果について、図５Ａ〜図６を参照して説明する。シミュレーション対象として撹拌槽モデルを仮定し、撹拌槽内の流体を撹拌するときの流体の挙動を解析した。
【００１３】
図５Ａは、シミュレーション対象の撹拌槽モデルの縦断面図であり、図５Ｂは図５Ａの一点鎖線５Ｂ−５Ｂにおける横断面図である。図５Ｂの一点鎖線５Ａ−５Ａにおける断面図が図５Ａに相当する。撹拌槽モデルは、撹拌槽１０及び撹拌槽１０に収容された流体１５を含む。
【００１４】
撹拌槽１０は、円筒状の側面の上端及び下端を２：１半楕円体の槽底で塞いだものである。槽径Ｄ１を０．３１ｍとした。撹拌槽１０内に、回転翼１１が配置されている。回転翼１１として、マックスブレンド（登録商標）に搭載される形状のものを用いた。翼径Ｄ２を０．１６ｍとした。撹拌槽１０の側面の内側に、周方向に等間隔で４枚の邪魔板１２を配置した。撹拌槽１０内に流体１５を収容した。
【００１５】
撹拌槽１０、回転翼１１、邪魔板１２、及び流体１５を、それぞれ粒子の集合体とし、この粒子系の挙動をシミュレーションした。撹拌槽１０、回転翼１１、及び邪魔板１２は剛体として取り扱った。撹拌槽１０を含む空間が解析領域となる。撹拌槽１０は剛体であるため、解析領域の境界に周期境界条件等を設定する必要はない。
【００１６】
シミュレーションで用いる重力加速度を実際の重力加速度より大きくし、重力加速度を変更してもフルード数が変化しないように回転翼１１の回転速度を変更した。実際のシミュレーションでは、回転翼１１の回転数を１３５０ｒｐｍとし、重力加速度を１４０６１１７．６５ｍ／ｓ２とした。回転翼１１の回転数が、フルード数を定義する特性速度に相当する。
【００１７】
図６は、従来の分子動力学法を用いて流体の挙動をシミュレーションした結果を示す図である。図６において、左から右に向かって時間が経過している。時間が経過するに従って流体の液面が下がり、流体の密度が高くなってしまっていることが確認される。液面が下がってしまうのは、重力加速度を実際の値より大きくしたためであると考えられる。
【００１８】
実施例によるシミュレーション方法では、流体の密度が高くなってしまわないように密度制御を行う。次に、図１を参照して、実施例によるシミュレーション方法について説明する。
【００１９】
図１は、実施例によるシミュレーション方法のフローチャートである。まず、シミュレーションの初期設定条件を決定する（ステップＳ１）。初期設定条件には、粒子系の各粒子の位置及び速度の初期値、各粒子の質量、粒子密度の設定値、設定温度、粒子の集合体で表された剛体の挙動（例えば回転翼の回転数等）、シミュレーションで用いる重力加速度、粒子間の相互作用を定義するパラメータ等が含まれる。
【００２０】
次に、粒子間の相互作用、すなわち各粒子に作用する力を計算する（ステップＳ２）。分子動力学計算で用いる粒子間の相互作用ポテンシャルとして、例えばレナードジョーンズポテンシャルを用いることができる。
【００２１】
相互作用を計算した後、流体を構成する複数の粒子の各々について、当該粒子を中心とする局所領域の粒子密度を算出する（ステップＳ３）。ここで、「粒子密度」は、単位体積当たりの流体の粒子の個数を意味する。粒子密度を算出するときの局所領域の大きさは、例えば、粒子間の相互作用の計算を行うときに計算を打ち切るカットオフ半径に基づいて決定するとよい。例えば、局所領域を、着目する粒子の位置を中心し、半径がカットオフ半径と等しい球体とするとよい。カットオフ半径として、例えば、粒子間の相互作用ポテンシャルの値が０になる距離の２．５倍〜５．５倍とするとよい。
【００２２】
粒子密度を算出した後、局所領域ごとに、粒子密度の算出値と設定値とに基づいて、設定温度を更新する（ステップＳ４）。例えば、粒子密度の算出値が設定値より高いとき、設定温度を現在の値より高くし、粒子密度の算出値が設定値より低いとき、設定温度を現在の値より低くする。設定温度は、粒子ごとに決定される。設定温度の変化量の決定方法については、後に詳しく説明する。
【００２３】
設定温度を更新した後、更新された設定温度で温度制御を行う（ステップＳ５）。例えば、局所領域の温度が設定温度に近づくように、粒子の運動エネルギを修正する。温度制御には、例えばランジュバン法を用いることができる。ランジュバン法においては、粒子にランダム力と散逸力とを与えて温度制御を行う。ランジュバン法については、Ｊ．Ｍ．テッセン著、「計算物理学」丸善出版（２００３年）の２０４〜２０８頁に説明されている。なお、温度制御に他の手法を適用してもよい。
【００２４】
温度制御を行った後、粒子の現在の状態に基づいて各粒子について運動方程式を数値的に解くことにより、粒子系の次の状態（１タイムステップ経過後の状態）を求める（ステップＳ６）。粒子系の状態は、複数の粒子の各々の位置及び速度で定義される。剛体を構成する複数の粒子については、粒子の相対位置関係が変化しないと仮定する。運動方程式を解くときに用いる各粒子に作用する力は、ステップＳ２で相互作用ポテンシャルに基づいて求められた力、及びステップＳ５の温度制御によって求められたランダム力及び散逸力を含む。この分子動力学計算により、１タイムステップ経過後の粒子の各々の位置及び速度が求まる。
【００２５】
粒子系の状態を更新した後、計算を終了するか否かを判定する（ステップＳ７）。例えば、計算を行ったステップ数が、予め決められた計算打ち切りステップ数に達したら、計算を終了する。または、粒子系の状態が定常状態に達したら計算を終了する。計算を終了しない場合は、ステップＳ２からステップＳ６までの処理を繰り返す。
【００２６】
次に、ステップＳ４において設定温度を更新する方法の一例について説明する。
粒子密度の設定値をρｃ、ｉ番目の粒子を中心とする局所領域の現在の粒子密度の算出値をρｉ、設定温度をＴｃで表したとき、以下の微分方程式を解くことによって、設定温度Ｔｃを変化させる。
【数１】
2019125296.tif 000003
ここで、ｍｉは粒子の質量、ｋＢはボルツマン定数、τは時定数、Ｔ０は温度の次元を持つパラメータ、Ｐ０は圧力の次元を持つパラメータである。
【００２７】
式（１）は、設定温度Ｔｃの加速度項と、粒子密度の算出値ρｉと設定値ρｃとの差に比例する外力と、粒子の質量ｍｉに比例する粘性摩擦力とを含む運動方程式と考えることができる。
【００２８】
次に、式（１）の導出過程について説明する。理想気体の状態方程式は以下の式で表される。
【数２】
2019125296.tif 000004
ここで、Ｐは圧力、Ｖは体積、ｎは物質量、Ｒは気体定数、Ｔは温度である。
【００２９】
粒子数をＮ、アボガドロ定数をＮＡで表すと、物質量ｎは、以下の式で表される。
【数３】
2019125296.tif 000005
【００３０】
式（２）及び（３）から、以下の式が導出される。
【数４】
2019125296.tif 000006
式（４）は、以下のように書き直すことができる。
【数５】
2019125296.tif 000007
ここで、ｍは粒子の質量、ρは単位体積当たりの粒子数（粒子密度）である。
【００３１】
式（５）から、粒子密度ρは以下の式で表される。
【数６】
2019125296.tif 000008
【００３２】
圧力Ｐが一定の条件の下で、式（６）から、以下の式が導出される。
【数７】
2019125296.tif 000009
ここで、ｄρは粒子密度変化量、ｄＴは温度の変化量である。式（７）を変形すると、式（８）が得られる。
【数８】
2019125296.tif 000010
【００３３】
式（７）の右辺のマイナス符号は、温度Ｔが高くなると粒子密度ρが低下し、温度Ｔが低下すると粒子密度ρが高くなることを意味する。
【００３４】
式（８）を、粒子密度の算出値ρｉが時定数τで設定値ρｃに近づくような運動方程式の形に書き換えると、以下の式が得られる。
【数９】
2019125296.tif 000011
粒子密度の算出値ρｉの振動に対しても算出値ρｉが設定値ρｃに近づくように減衰項を加えると、上記式（１）が得られる。
【００３５】
上記実施例によるシミュレーション方法を適用することにより得られる効果を確認するために、実際にシミュレーションを行った。図５Ａ及び図５Ｂに示した撹拌槽モデルをシミュレーション対象とした。シミュレーションにおいて、時間刻み幅ｄｔに対して時定数τを１００ｄｔ秒とし、温度の次元を持つパラメータＴ０を１０Ｋとし、圧力の次元を持つパラメータを２０００．０×１０５Ｐａとした。
【００３６】
次に、図２及び図３を参照して、撹拌槽モデルのシミュレーション結果について説明する。
図２は、シミュレーション結果を示す図である。図２において、左から右に向かって時間が経過している。濃い色の領域の粒子密度が、淡い色の領域の粒子密度より高い。時間が経過しても、流体の液面の高さが保たれていることがわかる。これは、流体の粒子密度が、ある程度設定値に維持されていることを意味する。
【００３７】
図３は、レイノルズ数Ｒｅと動力数Ｎｐとの関係を示すグラフである。横軸はレイノルズ数Ｒｅを表し、縦軸は動力数Ｎｐを表す。動力数Ｎｐは、回転翼１１（図５Ａ、図５Ｂ）に加わる動力を無次元化したものである。無次元化された動力数Ｎｐは、回転翼１１を構成する複数の粒子に作用する力に基づいて算出することができる。
【００３８】
図３のグラフにおいて、丸記号及び四角記号は、それぞれ実施例（図２）及び従来例（図６）によるシミュレーションにより得られたレイノルズ数Ｒｅ及び動力数Ｎｐを示す。加藤らによる「大型翼の撹拌所要動力に対する液高さの影響」（化学工学論文集、第４２巻、第６号、１８７頁〜１９１頁、２０１６年）に示された動力数とレイノルズ数との関係の予測式に従うグラフを実線で示す。
【００３９】
従来例による方法で得られた動力数Ｎｐは、予測式により予測された動力数Ｎｐから大幅に大きくなっている。従来例では、粒子密度が設定値より大きくなったことにより、流体の粘度が上昇し、回転翼１１にかかる動力が大きくなってしまっていると考えられる。これに対し、実施例による方法で得られた動力数Ｎｐと、予測式により予測された動力数Ｎｐとの誤差は５％程度であり、両者はよく一致している。実施例においては、重力加速度を実際の値より大きくしても、図１のステップＳ３〜Ｓ５の密度制御を行うことにより、十分な精度でシミュレーションを行うことが可能であることが確認された。
【００４０】
温度制御の設定温度を変化させることにより、粒子密度をほぼ設定値に保つことができる原理は、以下のように考えることができる。実施例によるシミュレーションにおいて、粒子密度が設定値より大きくなった領域においては、温度制御の設定温度が高くなる。流体の温度の上昇は、例えば式（７）に示したように、粒子密度を低下させる方向に作用する。逆に、粒子密度が設定値より小さくなった領域においては、温度制御の設定温度が低くなる。流体の温度の低下は、粒子密度を上昇させる方向に作用する。このように、粒子密度の算出値と設定値との差に基づいて温度制御の設定温度を更新することにより、粒子密度を設定値に維持するような作用が生じる。その結果、粒子密度をほぼ設定値に保ちつつシミュレーションを行うことができる。
【００４１】
上記実施例は、くりこみ群分子動力学、及び分子動力学のいずれにも適用することが可能である。また、上記実施例では粒子に作用する重力を実際の重力より大きくしたときに、フルード数が変化しないように特性速度を実際の特性速度より速くしたが、重力や特性速度等のシミュレーション条件を変化させる場合、必ずしもフルード数を不変にすることは必須ではない。シミュレーションの対象及びシミュレーションの目的によっては、シミュレーション条件を変化させる前後において、フルード数が変化してもよい。
【００４２】
次に、図４を参照して実施例によるシミュレーション装置について説明する。実施例によるシミュレーション装置は、図１に示したシミュレーション方法をコンピュータに実行させるプログラムを含み、このプログラムを実行する。
【００４３】
図４は、実施例によるシミュレーション装置のブロック図である。このシミュレーション装置として、シミュレーションプログラムを実行するコンピュータを用いることができる。このコンピュータは、中央処理装置（ＣＰＵ）２０、メインメモリ２１、補助記憶装置２２、入力装置２３、及び出力装置２４を含む。入力装置２３として、キーボード、ポインティングデバイス、通信装置（受信装置）、リムーバブルメディアの読み取り装置等の入力手段を用いることができる。出力装置２４として、ディスプレイ、プリンタ、通信装置（送信装置）、リムーバブルメディアへの書き込み装置等の出力手段を用いることができる。
【００４４】
補助記憶装置２２は、シミュレーションプログラムが格納されている記録媒体を含む。この記録媒体は、補助記憶装置２２に内蔵されたものでもよく、補助記憶装置２２に対して着脱可能なリムーバブルメディアでもよい。このシミュレーションプログラムがメインメモリ２１にロードされ、ＣＰＵ２０によって実行される。プログラム及びＣＰＵ２０を含む処理手段が、図１に示した各処理を実行する機能を持つ。プログラムは、コンピュータにこれらの機能を実現させる。
【００４５】
入力装置２３からシミュレーションに必要な情報が入力される。例えば、ステップＳ１（図１）で決定される種々の初期設定条件が入力される。なお、シミュレーション装置のユーザが入力する値は、これらの初期設定条件を構成するパラメータの値自体である必要はなく、これらの値を特定するためのパラメータでもよい。この場合、ＣＰＵ２０は、入力されたパラメータに基づいて、初期設定条件を算出する。
【００４６】
ＣＰＵ２０は、入力装置２３から入力された初期設定条件等を取得する。さらに、ＣＰＵ２０は、シミュレーション結果を出力装置２４に出力する。例えば、ステップＳ７（図１）で計算の終了と判定された場合、ＣＰＵ２０は、図２に示したような粒子の挙動を可視化した画像を出力装置２４に表示する。さらに、シミュレーションによって得られた動力数Ｎｐの値を出力装置２４に数値で表示する。
【００４７】
次に、上記実施例の変形例について説明する。上記実施例では、図１に示したように、粒子間の相互作用を計算（ステップＳ２）した後に、ステップＳ３からステップＳ５までの密度制御を行ったが、相互作用の計算を密度制御の後に実行してもよい。
【００４８】
上述の実施例は例示であり、本発明は上述の実施例に制限されるものではない。例えば、種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。
【符号の説明】
【００４９】
１０  撹拌槽
１１  回転翼
１２  邪魔板
１５  流体
２０  中央処理装置（ＣＰＵ）
２１  メインメモリ
２２  補助記憶装置
２３  入力装置
２４  出力装置
【図１】
2019125296.tif 000012
【図２】
2019125296.tif 000013
【図３】
2019125296.tif 000014
【図４】
2019125296.tif 000015
【図５】
2019125296.tif 000016
【図６】
2019125296.tif 000017