(57)【要約】【課題】高分子からなる流体の解析においても、壁面での流速がほぼゼロになり、実際の流速の分布を反映した解析を行うことが可能なシミュレーション装置を提供する。【解決手段】壁面に接触する流体を複数の流体粒子で表して、複数の流体粒子の挙動を解析する。壁面の形状を表す情報、壁面に配置する複数の壁面粒子を定義する情報、及び複数の流体粒子を定義する情報を取得する。取得された情報に基づいて、複数の壁面粒子と複数の流体粒子とが相互作用する条件の下で、複数の流体粒子の位置を時間発展させる。【選択図】図３
2020126531.tif 000002
【特許請求の範囲】
【請求項１】
壁面に接触する流体を複数の流体粒子で表して、前記複数の流体粒子の挙動を解析するシミュレーション装置であって、
前記壁面の形状を表す情報、前記壁面に配置する複数の壁面粒子を定義する情報、及び前記複数の流体粒子を定義する情報を取得し、
取得された情報に基づいて、前記複数の壁面粒子と前記複数の流体粒子とが相互作用する条件の下で、前記複数の流体粒子の位置を時間発展させるシミュレーション装置。
【請求項２】
前記流体は高分子材料を含み、前記複数の流体粒子のうち一部の流体粒子が結合して高分子を構成しており、前記壁面に沿って、前記流体粒子がすり抜けることができる隙間を持たせて１層目の壁面粒子を配置し、１層目より深い位置に、前記１層目の壁面粒子の隙間を塞ぐ位置に２層目の壁面粒子を配置する請求項１に記載のシミュレーション装置。
【請求項３】
壁面に接触する流体を複数の流体粒子で表して、前記複数の流体粒子の挙動を解析するシミュレーション方法であって、
前記壁面の形状を表す情報、前記壁面に配置する複数の壁面粒子を定義する情報、及び前記複数の流体粒子を定義する情報を決定し、
決定された情報に基づいて、前記複数の壁面粒子と前記複数の流体粒子とが相互作用する条件の下で、前記複数の流体粒子の位置を時間発展させるシミュレーション方法。
【請求項４】
前記流体は高分子材料を含み、前記複数の流体粒子のうち一部が結合して高分子を構成しており、
前記複数の流体粒子の位置を時間発展させる前に、前記壁面に沿って、前記流体粒子がすり抜けることができる隙間を持たせて１層目の壁面粒子を配置し、１層目より深い位置に、前記１層目の壁面粒子の隙間を塞ぐ位置に２層目の壁面粒子を配置する請求項３に記載のシミュレーション方法。
【請求項５】
壁面に接触する流体を複数の流体粒子で表して、前記複数の流体粒子の挙動を解析するシミュレーションをコンピュータで行うためのプログラムであって、
前記壁面の形状を表す情報、前記壁面に配置する複数の壁面粒子を定義する情報、及び前記複数の流体粒子を定義する情報を取得する機能と、
取得された情報に基づいて、前記複数の壁面粒子と前記複数の流体粒子とが相互作用する条件の下で、前記複数の流体粒子の位置を時間発展させる機能と
をコンピュータに実現させるためのプログラム。
【請求項６】
前記流体は高分子材料を含み、前記複数の流体粒子のうち一部が結合して高分子を構成しており、前記壁面に沿って、前記流体粒子がすり抜けることができる隙間を持たせて１層目の壁面粒子を配置し、１層目より深い位置に、前記１層目の壁面粒子の隙間を塞ぐ位置に２層目の壁面粒子を配置する機能を、さらにコンピュータに実現させるための請求項５に記載のプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、流体の流れを解析するシミュレーション装置、シミュレーション方法、及びプログラムに関する。
【背景技術】
【０００２】
壁面に接している流体の流れを、分子動力学法を用いて解析する手法が下記の非特許文献１に開示されている。壁面と接する位置（壁面境界）において流体の平均速度がゼロであることが妥当であることが実験から分かっている。非特許文献１に開示された手法では、壁面境界にミラー境界条件を適用し、壁面での流体の接線方向の平均速度がゼロになるように粒子の速度をリセットしている。
【先行技術文献】
【非特許文献】
【０００３】
【非特許文献１】T. Ishiwata, T. Murakami, S. Yukawa and N. Ito, “Particle Dynamics Simulations of the Navier-Stokes Flow with Hard Disks”, International Journal of Modern Physics C, Vol. 15, No.10 (2004), pp.1413-1424
【発明の概要】
【発明が解決しようとする課題】
【０００４】
流体を高分子とし、非特許文献１に記載されているミラー境界条件を適用して流体の流れを解析したところ、壁面境界で粒子の滑りが観測され、流速がゼロにならないことが判明した。なお、高分子の粒子の解析には、クレマーグレストモデル（Ｋｒｅｍｅｒ−Ｇｒｅｓｔモデル）を用いた。
【０００５】
本発明の目的は、高分子からなる流体の解析においても、壁面での流速がほぼゼロになり、実際の流速の分布を反映した解析を行うことが可能なシミュレーション装置、シミュレーション方法、及びプログラムを提供することである。
【課題を解決するための手段】
【０００６】
本発明の一観点によると、
壁面に接触する流体を複数の流体粒子で表して、前記複数の流体粒子の挙動を解析するシミュレーション装置であって、
前記壁面の形状を表す情報、前記壁面に配置する複数の壁面粒子を定義する情報、及び前記複数の流体粒子を定義する情報を取得し、
取得された情報に基づいて、前記複数の壁面粒子と前記複数の流体粒子とが相互作用する条件の下で、前記複数の流体粒子の位置を時間発展させるシミュレーション装置が提供される。
【０００７】
本発明の他の観点によると、
壁面に接触する流体を複数の流体粒子で表して、前記複数の流体粒子の挙動を解析するシミュレーション方法であって、
前記壁面の形状を表す情報、前記壁面に配置する複数の壁面粒子を定義する情報、及び前記複数の流体粒子を定義する情報を決定し、
決定された情報に基づいて、前記複数の壁面粒子と前記複数の流体粒子とが相互作用する条件の下で、前記複数の流体粒子の位置を時間発展させるシミュレーション方法が提供される。
【０００８】
本発明のさらに他の観点によると、
壁面に接触する流体を複数の流体粒子で表して、前記複数の流体粒子の挙動を解析するシミュレーションをコンピュータで行うためのプログラムであって、
前記壁面の形状を表す情報、前記壁面に配置する複数の壁面粒子を定義する情報、及び前記複数の流体粒子を定義する情報を取得する機能と、
取得された情報に基づいて、前記複数の壁面粒子と前記複数の流体粒子とが相互作用する条件の下で、前記複数の流体粒子の位置を時間発展させる機能と
をコンピュータに実現させるためのプログラムが提供される。
【発明の効果】
【０００９】
高分子からなる流体の解析においても、壁面での流速がほぼゼロになり、実際の流速の分布を反映した解析を行うことが可能である。
【図面の簡単な説明】
【００１０】
【図１】図１は、流体が流れる円管の中心軸を含む断面図である。
【図２】図２は、実施例によるシミュレーション装置のブロック図である。
【図３】図３Ａは、実施例によるシミュレーション装置でシミュレーション対象となる解析モデルの一例を示す図であり、図３Ｂは、円管の壁面及びその近傍に位置する高分子を示す模式図である。
【図４】図４は、円管の中心軸を含む断面において、壁面の近傍の流体粒子が感じるポテンシャルを示す図である。
【図５】図５は、実施例によるシミュレーション方法のフローチャートである。
【図６】図６Ａは、ニュートン流体の流れをシミュレーションしたときの流体粒子の速度分布のシミュレーション結果を示すグラフであり、図６Ｂは、べき乗則流体の流れをシミュレーションしたときの流体粒子の速度分布のシミュレーション結果を示すグラフである。
【発明を実施するための形態】
【００１１】
実施例について説明する前に、高分子からなる流体の円管内における流束分布について説明する。
【００１２】
図１は、流体が流れる円管１０の中心軸を含む断面図である。円管１０の中心軸に平行な方向をｚ軸方向とし、中心軸からの距離をｒで表す。円管１０の内側の半径をＲで表す。高分子材料においては、粘度とひずみ速度との関係がべき乗則に従うことが知られている。べき乗則流体のせん断応力をτ、粘性係数をη０、ひずみ速度をγドットで表すと、以下の式が成り立つ。
【数１】
2020126531.tif 000003
ここで、式（１）の右辺の
【数２】
2020126531.tif 000004
は、みかけの粘度となる。ｎは定数である。ニュートン流体においてはｎが１であり、高分子の流体においては、通常、ｎが１以下である。
【００１３】
円管内の流速分布には理論解が存在し、速度ｖ（ｒ）は以下の式で表される。
【数３】
2020126531.tif 000005
ここで、ρは流体の密度、ｇは重力加速度であり、ρｇは体積力を表す。速度ｖ（ｒ）は、円管の中心（ｒ＝０）で最大となり、円管の壁面（ｒ＝Ｒ）でゼロになる。
【００１４】
ニュートン流体の流れを、分子動力学法を用いてシミュレーションする際に、壁面境界にミラー境界条件を適用し、壁面での流体のｚ軸方向の速度がゼロになるように粒子の速度をリセットする処理を行うと、シミュレーションによって得られた速度分布は、壁面においてほぼゼロになる。ところが、式（１）の定数ｎが１以下のべき乗則流体の速度分布をシミュレーションによって求めると、壁面における速度がゼロにならないことが判明した。以下に説明する実施例では、高分子等のべき乗則流体のシミュレーションにおいても、壁面における粒子の速度がほぼゼロになる。
【００１５】
次に、図２〜図５Ｂを参照して、実施例によるシミュレーション装置及びシミュレーション方法について説明する。
【００１６】
図２は、実施例によるシミュレーション装置のブロック図である。実施例によるシミュレーション装置は、入力部３０、処理部３１、出力部３２、及び記憶部３３を含む。入力部３０から処理部３１にシミュレーション条件等が入力される。さらに、オペレータから入力部３０に各種指令（コマンド）等が入力される。入力部３０は、例えば通信装置、リムーバブルメディア読取装置、キーボード等で構成される。
【００１７】
処理部３１は、入力されたシミュレーション条件及び指令に基づいて分子動力学法またはくりこみ群分子動力学法（以下、単に分子動力学法という。）を用いたシミュレーションを行う。さらに、シミュレーション結果を出力部３２に出力する。シミュレーション結果には、シミュレーション対象物である粒子系の粒子の状態、粒子系の物理量の時間的変化等を表す情報が含まれる。処理部３１は、例えばコンピュータの中央処理ユニットを含む。分子動力学法によるシミュレーションをコンピュータに実行させるためのプログラムが、記憶部３３に記憶されている。出力部３２は、通信装置、リムーバブルメディア書込み装置、ディスプレイ等を含む。
【００１８】
図３Ａは、実施例によるシミュレーション装置で行うシミュレーションの解析モデルの一例を示す図である。円管１０の内側の壁面１１に接触する高分子材料からなる流体が、円管１０内を流れる。流体は、複数の流体粒子２１の集合体として表される。いくつかの流体粒子２１が結合して１つの高分子２０を構成している。流体は、複数の高分子２０を含む。流体粒子２１は、高分子２０を構成するモノマーに相当する。
【００１９】
図３Ｂは、円管１０の壁面１１及びその近傍に位置する高分子２０を示す模式図である。円管１０の壁面１１に沿って、複数の壁面粒子１５が配置される。壁面１１の近傍に位置する流体粒子２１は、壁面粒子１５の隙間に入り込む場合がある。
【００２０】
図３Ａ及び図３Ｂに示した解析モデルにおいて、高分子２０からなる粒子系の挙動を、分子動力学法により解析する。このとき、流体粒子２１の各々には、流体粒子２１同士の相互作用、及び流体粒子２１と壁面粒子１５との相互作用による力、及び外部から与えられる体積力が作用する。円管１０の軸方向両端の端面においては、周期境界条件を適用する。
【００２１】
次に、流体粒子２１同士の相互作用について説明する。
流体粒子２１同士の相互作用ポテンシャルφＦ（ｒ）として、任意の流体粒子２１の間においては、
【数４】
2020126531.tif 000006
が適用される。ここで、ｒは、流体粒子２１からの距離を表す。
【００２２】
ポテンシャルＵＦ０（ｒ）は、基本的に下記の式で表される。
【数５】
2020126531.tif 000007
ここで、ｆは無次元化された関数を表し、εＦ及びσＦは、流体粒子２１を特徴づけるフィッティングパラメータである。パラメータεＦはエネルギの次元を持ち、相互作用係数と呼ばれる。パラメータσＦは距離の次元を持ち、粒子の大きさに依存する。ポテンシャルＵＦ０（ｒ）として、例えばレナードジョーンズ型ポテンシャルを適用することができる。なお、その他に、モース型ポテンシャルを適用することも可能である。
【００２３】
同一の高分子２０内の相互に隣接する流体粒子２１の間においては、ポテンシャルＵＦ０（ｒ）に有限伸長非線形弾性ポテンシャルＵｃｈ（ｒ：εＦ，σＦ）が追加されて
【数６】
2020126531.tif 000008
が適用される。有限伸長非線形弾性ポテンシャルＵｃｈ（ｒ：εＦ，σＦ）は、ポテンシャルＵＦ０（ｒ）を定義するパラメータεＦ及びσＦに依存するパラメータを含む。
【００２４】
次に、流体粒子２１と壁面粒子１５との間の相互作用について説明する。
壁面粒子１５の間に作用する相互作用ポテンシャルφＷ（ｒ）として、
【数７】
2020126531.tif 000009
が適用される。
【００２５】
ポテンシャルＵＷ０（ｒ）は、基本的に下記の式で表される。
【数８】
2020126531.tif 000010
ここで、ｆは無次元化された関数を表し、εＷ及びσＷは、壁面粒子１５を特徴づけるフィッティングパラメータである。ポテンシャルＵＷ０（ｒ）として、例えばレナードジョーンズ型ポテンシャルを適用することができる。なお、その他に、モース型ポテンシャルを適用することも可能である。
【００２６】
流体粒子２１と壁面粒子１５との間の相互作用ポテンシャルφＦＷ（ｒ）として、ミキシングルールを適用して、
【数９】
2020126531.tif 000011
が適用される。ここで、ｒは流体粒子２１と壁面粒子１５との間の距離である。
【００２７】
次に、壁面粒子１５の配置方法について説明する。
まず、円管１０の三次元形状をテトラメッシュに分割する。内側の壁面１１の面上に位置する節点の位置に、１層目の壁面粒子１５を配置する。このとき、１層目の壁面粒子１５の隙間を、流体粒子２１がすり抜けることができる程度の大きさとする。
【００２８】
１層目の壁面粒子１５を配置した後、１層目よりも深い位置に２層目の複数の壁面粒子１５を配置する。２層目の壁面粒子１５は、１層目の壁面粒子１５の隙間を埋める位置に配置する。例えば、１層目の、相互に隣り合う３つの壁面粒子１５の位置を３つの頂点とするほぼ正四面体の残りの１つの頂点の位置に、２層目の壁面粒子１５を配置する。２層目の壁面粒子１５が、１層目の壁面粒子１５の隙間を埋めているため、流体粒子２１は、１層目及び２層目の壁面粒子１５からなる壁をすり抜けることはできない。
【００２９】
図４は、円管１０の中心軸を含む断面において、壁面１１の近傍の流体粒子２１が感じるポテンシャルの計算結果を示す図である。図４に示した破線は、等ポテンシャル面１６を示す。壁面粒子１５により、凹凸を持つ等ポテンシャル面１６が形成されていることがわかる。
【００３０】
図５は、実施例によるシミュレーション方法のフローチャートである。図５に示した各処理は、処理部３１（図２）が記憶部３３（図２）に格納されているプログラムを実行することにより実現される。
【００３１】
まず、シミュレーション対象の円管１０の壁面１１の形状、壁面粒子１５及び流体粒子２１を定義する情報、初期条件、その他のシミュレーション条件を決定する。これらの情報を入力部３０に入力する。処理部３１は、入力部３０から、シミュレーション対象の円管１０の壁面１１の形状を定義する形状定義データ、壁面粒子１５及び流体粒子２１を定義する情報、初期条件、及びその他の必要なシミュレーション条件を取得する（ステップＳ１）。
【００３２】
壁面粒子１５及び流体粒子２１を定義する情報には、例えば、式（５）及び式（８）のフィッティングパラメータの値、式（６）の有限伸長非線形弾性ポテンシャルＵｃｈ（ｒ：εＦ，σＦ）を定義する情報、粒子の質量等が含まれる。初期条件には、流体粒子２１の位置や速度の初期値を定義する情報が含まれる。その他のシミュレーション条件には、流体に作用する体積力を規定するための流体の密度や重力に関する情報、流体の粘性係数に関する情報等が含まれる。
【００３３】
壁面１１の形状に基づいて、壁面１１に沿って複数の壁面粒子１５を配置する（ステップＳ２）。壁面粒子１５の配置には、図３を参照して説明した手法を用いる。次に、円管１０の内側の空間に流体粒子２１を配置する（ステップＳ３）。
【００３４】
壁面粒子１５と流体粒子２１との相互作用、流体粒子２１同士の相互作用、及び流体粒子２１に作用する体積力を考慮して運動方程式を解くことにより、流体粒子２１の位置を時間発展させる（ステップＳ４）。このとき、壁面粒子１５の位置は不変とする。なお、壁面粒子１５の間に有限伸長非線形弾性ポテンシャルによる力が作用していると仮定して、壁面粒子１５についても運動方程式を解き、その位置を時間発展させてもよい。
【００３５】
予め決められた所定のタイムステップ数だけ、または流体粒子２１が定常状態に達するまで、ステップＳ４の計算を繰り返す（ステップＳ５）。
【００３６】
解析が終了すると、解析結果を出力部３２に出力する（ステップＳ６）。出力される情報には、例えば流体粒子２１の速度分布を表す情報等が含まれる。
【００３７】
次に、図６Ａ及び図６Ｂを参照して、実施例によるシミュレーション方法、及び壁面にミラー境界条件を適用した従来の手法を用いて流体の流れをシミュレーションした結果について説明する。
【００３８】
シミュレーション対象の円管１０（図３）の直径を４０ｍｍ、長さを８０ｍｍとした。長さ方向の両端面には周期境界条件を適用した。また、流体粒子２１に軸方向の体積力を加え、平均流速が２ｍ／ｓになるまで解析を行った。高分子２０の挙動の解析には、クレマーグレストモデルを適用した。
【００３９】
図６Ａは、ニュートン流体の流れをシミュレーションしたときの流体粒子の速度分布のシミュレーション結果を示すグラフである。横軸は円管１０の中心軸から半径方向への距離を単位「ｍｍ」で表し、縦軸は流体粒子２１のｚ方向（軸方向）の速度を単位「ｍ／ｓ」で表す。図６Ａに示したグラフにおいて、ハッチング付きの丸記号は、実施例によるシミュレーション方法を用いて解析した結果を示し、中空の丸記号は、ミラー境界条件を適用して解析した結果を示す。グラフ中の実線は、式（３）から求まる理論解を示す。
【００４０】
流体がニュートン流体である場合には、実施例による方法及びミラー境界条件を適用した方法のいずれにおいても、ほぼ理論解に近い速度分布が得られていることがわかる。
【００４１】
図６Ｂは、べき乗則流体の流れをシミュレーションしたときの流体粒子の速度分布のシミュレーション結果を示すグラフである。図６Ａの場合と同様に、図６Ｂ横軸は円管１０の中心軸から半径方向への距離を単位「ｍｍ」で表し、縦軸は流体粒子２１のｚ方向（軸方向）の速度を単位「ｍ／ｓ」で表す。図６Ｂに示したグラフにおいて、ハッチング付きの丸記号は、実施例によるシミュレーション方法を用いて解析した結果を示し、中空の丸記号は、ミラー境界条件を適用して解析した結果を示す。グラフ中の実線は、式（３）から求まる理論解を示す。
【００４２】
ミラー境界条件を適用した場合は、壁面からごくわずか離れた位置で、流体粒子２１の速度が約２ｍ／ｓになっている。ミラー境界条件を適用し、かつ壁面における軸方向の速度がゼロになるように速度をリセットして滑りなし条件を実現しているにも関わらず、壁面の極近傍において滑りが発生していることがわかる。また、流体粒子２１の速度は、径方向の位置に依らず、ほぼ一定であることがわかる。これに対し、実施例によるシミュレーション方法を用いた場合には、理論解とほぼ近いシミュレーション結果が得られており、壁面の極近傍において滑りが抑制されていることがわかる。
【００４３】
以下、上述のシミュレーション結果について考察する。
べき乗則流体においては、壁面１１において接線方向の速度がゼロになるようにリセットしたとしても、時間発展の次のタイムステップにおいて、壁面１１の極近傍の流体粒子２１が円管１０の中心軸側の流れに引きずられて軸方向に移動してしまう。このため、壁面１１の極近傍で滑りが生じてしまう。これに対し、実施例においては、壁面１１に等ポテンシャル面の凹凸が形成され、その凹部に流体粒子２１がトラップされる。凹部にトラップされた流体粒子２１は、時間発展の次のタイムステップにおいても軸方向に移動しにくいため、壁面１１の極近傍において滑りを抑制する効果が得られる。
【００４４】
次に、図６Ａ及び図６Ｂに示したシミュレーション結果からわかる本実施例の優れた効果について説明する。
【００４５】
ミラー境界条件を適用し、壁面における滑りなし条件を実現する方法は、ニュートン流体のシミュレーションに適用することが可能であるが、べき乗則流体に適用することはできない。これに対し、実施例によるシミュレーション方法は、ニュートン流体及びべき乗則流体の両方に適用することが可能である。
【００４６】
さらに、ミラー境界条件を適用する場合は、壁面において流体粒子の速度をリセットするため、流体粒子２１の全エネルギが保存されない。これに対し、実施例によるシミュレーション方法を適用する場合は、速度のリセットを行わないため、流体粒子２１の全エネルギが保存されるという優れた効果が得られる。
【００４７】
上記実施例では、円管の中を流れる流体の流れをシミュレーションしているが、実施例によるシミュレーション装置及びシミュレーション方法は、その他の形状の壁面に接する流体の流れのシミュレーションに適用することも可能である。例えば、実施例によるシミュレーション方法は、高分子材料を混錬機で混ぜるときの高分子材料の挙動のシミュレーションに適用することができる。
【００４８】
上述の実施例は例示であり、本発明は上述の実施例に制限されるものではない。例えば、種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。
【符号の説明】
【００４９】
１０  円管
１１  壁面
１５  壁面粒子
１６  等ポテンシャル面
２０  高分子
２１  流体粒子
３０  入力部
３１  処理部
３２  出力部
３３  記憶部
【図１】
2020126531.tif 000012
【図２】
2020126531.tif 000013
【図３】
2020126531.tif 000014
【図４】
2020126531.tif 000015
【図５】
2020126531.tif 000016
【図６】
2020126531.tif 000017