JP7285223B2 - シミュレーション装置、シミュレーション方法、及びプログラム - Google Patents

Description

本発明は、流体の流れを解析するシミュレーション装置、シミュレーション方法、及びプログラムに関する。

壁面に接している流体の流れを、分子動力学法を用いて解析する手法が下記の非特許文献１に開示されている。壁面と接する位置（壁面境界）において流体の平均速度がゼロであることが妥当であることが実験から分かっている。非特許文献１に開示された手法では、壁面境界にミラー境界条件を適用し、壁面での流体の接線方向の平均速度がゼロになるように粒子の速度をリセットしている。

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

流体を高分子とし、非特許文献１に記載されているミラー境界条件を適用して流体の流れを解析したところ、壁面境界で粒子の滑りが観測され、流速がゼロにならないことが判明した。なお、高分子の粒子の解析には、クレマーグレストモデル（Ｋｒｅｍｅｒ－Ｇｒｅｓｔモデル）を用いた。

本発明の目的は、高分子からなる流体等の解析においても、壁面での流速がほぼゼロになり、実際の流速の分布を反映した解析を行うことが可能なシミュレーション装置、シミュレーション方法、及びプログラムを提供することである。

本発明の一観点によると、
壁面に接触する流体を複数の流体粒子で表した解析モデルにおいて、前記複数の流体粒子の挙動を解析するシミュレーション装置であって、
前記壁面の形状を定義する情報、前記複数の流体粒子が前記壁面から受ける相互作用ポテンシャルを定義する情報、及び前記流体の物性値を含むシミュレーション条件が入力される入力部と、
前記入力部に入力されたシミュレーション条件を取得し、取得された情報に基づいて前記複数の流体粒子について運動方程式を解いて前記複数の流体粒子の位置を時間発展させる処理部と
を有し、
前記処理部は、
前記複数の流体粒子のうち前記壁面までの距離が近接判定閾値以下になった流体粒子を壁面近接粒子として検出し、
前記壁面近接粒子と相互作用を及ぼし合う位置に複数の仮想粒子を発生させ、前記複数の仮想粒子の位置は固定し、前記壁面近接粒子と前記複数の仮想粒子との間に、前記壁面に対して平行な方向への前記壁面近接粒子の移動を妨げる相互作用ポテンシャルを作用させて、前記壁面近接粒子について運動方程式を解くシミュレーション装置が提供される。

本発明の他の観点によると、
壁面に接触する流体を複数の流体粒子で表した解析モデルにおいて、前記複数の流体粒子の挙動を解析するシミュレーション方法であって、
前記壁面の形状を定義する情報、前記複数の流体粒子が前記壁面から受ける相互作用ポテンシャルを定義する情報、及び前記流体の物性値を含むシミュレーション条件を取得し、
取得された情報に基づいて前記複数の流体粒子について運動方程式を解いて前記複数の流体粒子の挙動を解析し、
解析中に、前記複数の流体粒子のうち前記壁面までの距離が近接判定閾値以下になった流体粒子を壁面近接粒子として検出し、
検出された前記壁面近接粒子と相互作用を及ぼし合う位置に複数の仮想粒子を発生させ、前記複数の仮想粒子の位置は固定し、前記壁面近接粒子と前記複数の仮想粒子の間に、前記壁面近接粒子の前記壁面に対して平行な方向への移動を妨げる相互作用ポテンシャルを作用させて、前記壁面近接粒子について運動方程式を解くシミュレーション方法が提供される。

本発明のさらに他の観点によると、
壁面に接触する流体を複数の流体粒子で表した解析モデルにおいて、前記複数の流体粒子の挙動を解析するシミュレーションをコンピュータに実行させるプログラムであって、
前記壁面の形状を定義する情報、前記複数の流体粒子が前記壁面から受ける相互作用ポテンシャルを定義する情報、及び前記流体の物性値を含むシミュレーション条件を取得する機能と、
取得された情報に基づいて前記複数の流体粒子について運動方程式を解いて前記複数の流体粒子の挙動を解析する機能と、
解析中に、前記複数の流体粒子のうち前記壁面までの距離が近接判定閾値以下になった流体粒子を壁面近接粒子として検出する機能と、
検出された前記壁面近接粒子と相互作用を及ぼし合う位置に複数の仮想粒子を発生させ、前記複数の仮想粒子の位置は固定し、前記壁面近接粒子と前記複数の仮想粒子の間に、前記壁面近接粒子の前記壁面に対して平行な方向への移動を妨げる相互作用ポテンシャルを作用させて、前記壁面近接粒子について運動方程式を解く機能と
をコンピュータに実現させるプログラムが提供される。

高分子からなる流体等の解析においても、壁面での流速がほぼゼロになり、実際の流速の分布を反映した解析を行うことが可能である。

図１は、流体が流れる円管の中心軸を含む断面図である。 図２は、実施例によるシミュレーション装置のブロック図である。 図３Ａは、実施例によるシミュレーション装置で行うシミュレーションの解析モデルの一例を示す図であり、図３Ｂは、流体粒子の間で作用する相互作用ポテンシャルφ_Ｆ（ｒ）の一例を示すグラフであり、図３Ｃは、壁面が流体粒子に及ぼす相互作用ポテンシャルφ_Ｗ（ｒ）の一例を示すグラフである。 図４は、本実施例によるシミュレーション方法のフローチャートである。 図５は、符号付距離関数を説明するための模式図である。 図６は、ステップＳ４（図４）の詳細な処理を示すフローチャートである。 図７は、壁面と流体粒子との位置関係を示す模式図である。 図８Ａは、壁面近接粒子から壁面下した垂線に対して平行な方向から見たときの壁面近接粒子及び複数の仮想粒子の位置関係を示す図であり、図８Ｂは、壁面近接粒子から壁面に下した垂線に対して垂直な方向から見たときの壁面近接粒子及び複数の仮想粒子の位置関係を示す図であり、図８Ｃは、相互作用ポテンシャルφ_ｖを示すグラフである。 図９は、べき乗則流体の流れをシミュレーションしたときのシミュレーション結果を示すグラフである。

実施例について説明する前に、高分子からなる流体の円管内における流速分布について説明する。

図１は、流体が流れる円管１０の中心軸を含む断面図である。円管１０の中心軸に平行な方向をｚ軸方向とし、中心軸からの距離をｒで表す。円管１０の内側の半径をＲで表す。高分子材料においては、粘度とひずみ速度との関係がべき乗則に従うことが知られている。べき乗則流体のせん断応力をτ、粘性係数をη_０、ひずみ速度をγドットで表すと、以下の式が成り立つ。

ここで、式（１）の右辺の

は、みかけの粘度となる。ｎは定数である。ニュートン流体においてはｎが１であり、高分子の流体においては、通常、ｎが１以下である。

円管内の流速分布には理論解が存在し、速度ｖ（ｒ）は以下の式で表される。

ここで、ρは流体の密度、ｇは重力加速度であり、ρｇは体積力を表す。速度ｖ（ｒ）は、円管の中心（ｒ＝０）で最大となり、円管の壁面（ｒ＝Ｒ）でゼロになる。

ニュートン流体の流れを、分子動力学法を用いてシミュレーションする際に、壁面境界にミラー境界条件を適用し、壁面での流体のｚ軸方向の速度がゼロになるように粒子の速度をリセットする処理を行うと、シミュレーションによって得られた速度分布は、壁面においてほぼゼロになる。ところが、式（１）の定数ｎが１以下のべき乗則流体の速度分布をシミュレーションによって求めると、壁面における速度がゼロにならないことが判明した。以下に説明する実施例では、高分子等のべき乗則流体のシミュレーションにおいても、壁面における粒子の速度がほぼゼロになる。

次に、図２～図８Ｃを参照して、実施例によるシミュレーション装置及びシミュレーション方法について説明する。

図２は、実施例によるシミュレーション装置のブロック図である。実施例によるシミュレーション装置は、入力部３０、処理部３１、出力部３２、及び記憶部３３を含む。入力部３０から処理部３１にシミュレーション条件等が入力される。さらに、オペレータから入力部３０に各種指令（コマンド）等が入力される。入力部３０は、例えば通信装置、リムーバブルメディア読取装置、キーボード等で構成される。

処理部３１は、入力されたシミュレーション条件及び指令に基づいて分子動力学法またはくりこみ群分子動力学法（以下、単に分子動力学法という。）を用いたシミュレーションを行う。さらに、シミュレーション結果を出力部３２に出力する。シミュレーション結果には、シミュレーション対象物である粒子系の粒子の状態、粒子系の物理量の時間的変化等を表す情報が含まれる。処理部３１は、例えばコンピュータの中央処理ユニット（ＣＰＵ）を含む。分子動力学法によるシミュレーションをコンピュータに実行させるためのプログラムが、記憶部３３に記憶されている。出力部３２は、通信装置、リムーバブルメディア書込み装置、ディスプレイ等を含む。

図３Ａは、実施例によるシミュレーション装置で行うシミュレーションの解析モデルの一例を示す図である。円管１０の内側の壁面１１に接触する高分子材料からなる流体が、円管１０内を流れる。流体は、複数の流体粒子２１の集合体として表される。いくつかの流体粒子２１が結合して１つの高分子２０を構成している。流体は、複数の高分子２０を含む。流体粒子２１は、高分子２０を構成するモノマーに相当する。

図３Ａに示した解析モデルにおいて、高分子２０からなる粒子系の挙動を、分子動力学法により解析する。このとき、流体粒子２１の各々には、流体粒子２１同士の相互作用による力、流体粒子２１と壁面１１との相互作用による力、及び外部から与えられる体積力が作用する。円管１０の軸方向に関して両端の端面においては、周期境界条件を適用する。

次に、流体粒子２１同士の相互作用について説明する。
流体粒子２１同士の相互作用ポテンシャルφ_Ｆ（ｒ）として、任意の流体粒子２１の間においては、

が適用される。ここで、ｒは、流体粒子２１からの距離を表す。

ポテンシャルＵ_Ｆ０（ｒ）は、基本的に下記の式で表される。

ここで、ｆは無次元化された関数を表し、ε_Ｆ及びσ_Ｆは、流体粒子２１を特徴づけるフィッティングパラメータである。フィッティングパラメータε_Ｆはエネルギの次元を持ち、相互作用係数と呼ばれる。フィッティングパラメータσ_Ｆは距離の次元を持ち、粒子の大きさに依存する。ポテンシャルＵ_Ｆ０（ｒ）として、例えばレナードジョーンズ型ポテンシャルを適用することができる。なお、その他に、モース型ポテンシャルを適用することも可能である。

図３Ｂは、流体粒子２１の間で作用する相互作用ポテンシャルφ_Ｆ（ｒ）の一例を示すグラフである。流体粒子２１の近傍において距離ｒが大きくなるにしたがってポテンシャルが低下し、流体粒子２１からの距離ｒがｒ_０の位置でポテンシャルが最小値を示す。流体粒子２１からの距離ｒがｒ_０より遠い範囲では、距離ｒが大きくなるにしたがってポテンシャルが徐々に大きくなり、０に漸近する。

同一の高分子２０内の相互に隣接する流体粒子２１の間においては、ポテンシャルＵ_Ｆ０（ｒ）に有限伸長非線形弾性ポテンシャルＵ_ｃｈ（ｒ：ε_Ｆ，σ_Ｆ）が追加されて

が適用される。有限伸長非線形弾性ポテンシャルＵ_ｃｈ（ｒ：ε_Ｆ，σ_Ｆ）は、ポテンシャルＵ_Ｆ０（ｒ）を定義するフィッティングパラメータε_Ｆ及びσ_Ｆに依存するパラメータを含む。

次に、流体粒子２１と壁面１１との間の相互作用について説明する。
流体粒子２１が壁面１１から受ける相互作用ポテンシャルφ_Ｗ（ｒ）として、

を適用する。

ポテンシャルＵ_Ｗ０（ｒ）は、ポテンシャルＵ_Ｆ０（ｒ）と同様に、基本的に下記の式で表される。

ここで、ｆは無次元化された関数を表し、ε_Ｗ及びσ_Ｗは、壁面１１を特徴づけるフィッティングパラメータである。ポテンシャルＵ_Ｗ０（ｒ）として、例えばレナードジョーンズ型ポテンシャルを適用することができる。なお、その他に、壁面１１に近づくにしたがって流体粒子２１に加わる斥力が大きくなるようなポテンシャル、例えばモース型ポテンシャル等を適用することも可能である。

図３Ｃは、壁面１１が流体粒子２１に及ぼす相互作用ポテンシャルφ_Ｗ（ｒ）の一例を示すグラフである。相互作用ポテンシャルφ_Ｗ（ｒ）の形状は、相互作用ポテンシャルφ_Ｆ（ｒ）の形状と類似している。本実施例では、相互作用ポテンシャルφ_Ｗ（ｒ）と、相互作用ポテンシャルφ_Ｆ（ｒ）とを同一としている。すなわち、式（８）のフィッティングパラメータε_Ｗ及びσ_Ｗが、それぞれ式（５）のフィッティングパラメータε_Ｆ及びσ_Ｆに等しい。なお、式（８）のフィッティングパラメータε_Ｗ及びσ_Ｗを、それぞれ式（５）のフィッティングパラメータε_Ｆ及びσ_Ｆに等しくする必要はない。フィッティングパラメータε_Ｗ及びσ_Ｗが異なる種々の条件でシミュレーションを行い、現実の流速分布を最もよく反映している条件を、フィッティングパラメータε_Ｗ及びσ_Ｗの値として採用してもよい。

図４は、本実施例によるシミュレーション方法のフローチャートである。図４に示した各処理は、処理部３１（図２）が記憶部３３（図２）に格納されているプログラムを実行することにより実現される。

まず、シミュレーション対象の円管１０の壁面１１（図３Ａ）の形状、流体粒子２１を定義する情報（例えば流体の物性値）、初期条件、その他のシミュレーション条件を決定する。これらの情報を入力部３０に入力する。処理部３１は、入力部３０から、シミュレーション対象の円管１０の壁面１１の形状を定義する形状定義データ、流体粒子２１を定義する情報、初期条件、及びその他の必要なシミュレーション条件を取得する（ステップＳ１）。

流体粒子２１を定義する情報には、例えば、式（５）のフィッティングパラメータε_Ｆ及びσ_Ｆの値、式（６）の有限伸長非線形弾性ポテンシャルＵ_ｃｈ（ｒ：ε_Ｆ，σ_Ｆ）を定義する情報、粒子の質量等が含まれる。本実施例においては、式（８）のフィッティングパラメータε_Ｗ及びσ_Ｗの値は、式（５）のフィッティングパラメータε_Ｆ及びσ_Ｆの値と同一とする。

初期条件には、流体粒子２１の位置や速度の初期値を定義する情報が含まれる。その他のシミュレーション条件には、流体に作用する体積力を規定するための流体の密度や重力に関する情報、流体の粘性係数に関する情報等が含まれる。

次に、壁面１１の形状に基づいて、符号付距離関数（ＳＤＦ）を生成する（ステップＳ２）。図５を参照して、符号付距離関数について説明する。

図５は、符号付距離関数を説明するための模式図である。流体の流路を含む空間を三次元直交格子で分割する。複数の格子点ＧＰの各々について、格子点ＧＰから壁面１１に下した垂線ＰＬの長さｒ_Ｗ（すなわち格子点ＧＰから壁面１１までの距離）を、格子点ＧＰの各々に対応させる。壁面１１より内側の格子点ＧＰに対しては正の距離を対応付け、外側の格子点ＧＰに対しては負の距離を対応付ける。図５では、空間を二次元で表しているが、実際には三次元の空間内の各格子点ＧＰに、壁面１１までの符号付の距離が対応付けられる。空間内の各位置に、壁面１１までの符号付の距離を対応付けた情報は、符号付距離関数といわれる。

符号付距離関数を用い、必要に応じて補間演算を行うことにより、空間内の任意の点について、壁面１１までの距離、及び壁面１１に下した垂線の方向を求めることができる。任意の点について、壁面１１までの距離と垂線の方向とがわかれば、壁面１１から受ける相互作用ポテンシャルに基づいて、流体粒子２１が壁面１１から受ける力を計算することができる。

図４のステップＳ２において符号付距離関数を生成した後、円管１０（図３Ａ）の壁面１１の内側の空間に、複数の高分子２０を構成する複数の流体粒子２１を配置する（ステップＳ３）。複数の流体粒子２１を配置した後、流体粒子２１ごとに、流体粒子２１に作用する相互作用ポテンシャルに基づいて運動方程式を解くことにより、流体粒子２１の位置を時間発展させる（ステップＳ４）。ステップＳ４の詳細な処理については、後に図６～図８Ｃを参照して説明する。解析終了条件を満たすまで、流体粒子２１の位置を時間発展させる処理を繰り返す（ステップＳ５）。例えば、流れ場が定常状態に達したら、解析終了条件が満たされたと判定する。解析が終了すると、解析結果を出力部３２に出力する（ステップＳ６）。出力される情報には、例えば流体粒子２１の速度分布を表す情報等が含まれる。

図６は、ステップＳ４（図４）の詳細な処理を示すフローチャートである。図６のステップＳ４１１からステップＳ４１９までの処理を、全ての流体粒子２１のそれぞれについて実行する（ステップＳ４２０）。

まず、着目する流体粒子２１が壁面１１に近接しているか否かを判定する（ステップＳ４１１）。

図７を参照して、ステップＳ４１１の処理について説明する。
図７は、壁面１１と流体粒子２１との位置関係を示す模式図である。流体粒子２１から壁面１１までの距離が近接判定閾値Ｌ１以下である場合、その流体粒子２１は壁面１１に近接していると判定する。本明細書において、壁面１１までの近接判定閾値Ｌ１以下の流体粒子２１を「壁面近接粒子」ということとする。

近接判定閾値Ｌ１として、例えば、式（７）及び図３Ｃに示した壁面１１と流体粒子２１との相互作用ポテンシャルφ_Ｗが最小となるときの距離ｒ_０を採用する。これは、流体粒子２１が壁面１１から斥力を受け始めたら、その流体粒子２１を壁面近接粒子２１Ａとして取り扱うことを意味する。物理的には、判定対象の流体粒子２１と壁面１１との間に他の流体粒子２１が存在しない場合、その流体粒子２１を壁面近接粒子２１Ａとして取り扱うことに相当する。

着目する流体粒子２１が壁面近接粒子であるとき、着目する流体粒子２１の位置を時間発展させる前（直近のタイムステップ実行前）においても、着目する流体粒子２１が壁面１１に近接していたか否かを判定する（ステップＳ４１２）。時間発展前には、着目する流体粒子２１が壁面１１に近接していなかった場合、すなわち、直近のタイムステップにおける流体粒子２１の移動により、壁面１１に近接していない位置から近接した位置に移動した場合には、着目する流体粒子２１の近傍に複数の仮想粒子を配置する（ステップＳ４１５）。

図８Ａ～図８Ｃを参照して、仮想粒子を配置する処理について説明する。
図８Ａは、壁面近接粒子２１Ａから壁面１１下した垂線に対して平行な方向から見たときの壁面近接粒子２１Ａ及び複数の仮想粒子２５の位置関係を示す図である。図８Ｂは、壁面近接粒子２１Ａから壁面１１に下した垂線２６に対して垂直な方向から見たときの壁面近接粒子２１Ａ及び複数の仮想粒子２５の位置関係を示す図である。

例えば、壁面近接粒子２１Ａから壁面１１に下した垂線２６に対して直交し、壁面近接粒子２１Ａを通過する仮想的な平面２７の上に、３個の仮想粒子２５を配置する。３個の仮想粒子２５は、それぞれ壁面近接粒子２１Ａの位置を重心とする正三角形の３つの頂点の位置に配置する。仮想的な平面２７の面内の回転方向に関する正三角形の姿勢はランダムとする。

次に、仮想粒子２５が壁面近接粒子２１Ａに及ぼす相互作用ポテンシャルφ_ｖについて説明する。相互作用ポテンシャルφ_ｖを、以下のように定義する。

ここで、ｒは壁面近接粒子２１Ａと仮想粒子２５との距離、ε_ｖ、σ_ｖはフィッティングパラメータである。式（９）の１行目の式の右辺の定数０．２５をゼロに置き換えると、相互作用ポテンシャルφ_ｖはレナードジョーンズ型のポテンシャルになる。

図８Ｃは、相互作用ポテンシャルφ_ｖを示すグラフである。壁面近接粒子２１Ａと仮想粒子２５との距離ｒがＬ_ｒｍ未満のとき、壁面近接粒子２１Ａが仮想粒子２５から斥力を受ける。距離ｒがＬ_ｒｍ以上のとき、壁面近接粒子２１Ａは仮想粒子２５から力を受けない。ここで、Ｌ_ｒｍは、式（９）の２^１／６σ_ｖに等しい。本明細書において、距離Ｌ_ｒｍを「斥力発生最長距離」ということとする。

仮想粒子２５（図８Ａ）を配置するとき、仮想粒子２５の各々と壁面近接粒子２１Ａとの距離は、斥力発生最長距離Ｌ_ｒｍと等しくする。このとき、壁面近接粒子２１Ａの位置において、仮想粒子２５による相互作用ポテンシャルφ_ｖの大きさはゼロである。従って、仮想粒子２５を配置する前後において系全体のエネルギは変化しない。

３個の仮想粒子２５の位置を頂点とする正三角形の内部では、正三角形の重心位置において相互作用ポテンシャルφ_ｖが最小になる。すなわち、壁面近接粒子２１Ａが正三角形の内部で移動すると、壁面近接粒子２１Ａを重心位置に押し戻そうとする力が壁面近接粒子２１Ａに作用する。言い換えると、３個の仮想粒子２５による相互作用ポテンシャルφ_ｖは、壁面近接粒子２１Ａの、壁面１１に対して平行な方向への移動を妨げる方向に作用する。なお、壁面近接粒子２１Ａが１タイムステップの時間発展によって正三角形の外側まで移動しないように、時間刻み幅を小さく設定することが好ましい。

相互作用ポテンシャルφ_ｖとして、式（９）のポテンシャル以外のものを採用してもよい。例えば、流体粒子２１から仮想粒子２５までの距離が斥力発生最長距離Ｌ_ｒｍ未満のとき、流体粒子２１に斥力が発生し、斥力発生最長距離Ｌ_ｒｍ以上のときには流体粒子２１が仮想粒子２５から力を受けないようなポテンシャルを採用してもよい。

図６のステップＳ４１５で仮想粒子２５を配置した後、壁面近接粒子２１Ａに作用する力を計算する（ステップＳ４１６）。壁面近接粒子２１Ａには、仮想粒子２５（図８Ａ、図８Ｂ）からの相互作用ポテンシャルφ_ｖ（図８Ｃ）、壁面からの相互作用ポテンシャルφ_Ｗ（図３Ｃ）、及び他の流体粒子２１からの相互作用ポテンシャルφ_Ｆ（図３Ｂ）による力が作用する。その後、着目する流体粒子２１について運動方程式を解き、その位置を時間発展させる（ステップＳ４１９）。

ステップＳ４１２において、時間発展前にも、着目する流体粒子２１が壁面１１に近接していたと判定された場合、着目する流体粒子２１の近傍には、既に３個の仮想粒子２５（図８Ａ、図８Ｂ）が配置されている。この仮想粒子２５の位置を固定した条件の下で、着目する流体粒子２１に作用する力を計算する（ステップＳ４１６）。その後、着目する流体粒子２１について運動方程式を解き、その位置を時間発展させる（ステップＳ４１９）。

ステップＳ４１１において、着目する流体粒子２１が壁面１１に近接していないと判定された場合、着目する流体粒子２１に仮想粒子２５が対応付けられているか否かを判定する（ステップＳ４１３）。着目する流体粒子２１が、仮想粒子２５に対応付けられていない場合、相互作用ポテンシャルφ_Ｗ（図３Ｃ）及びφ_Ｆに（図３Ｂ）基づいて、流体粒子２１に作用する力を計算する（ステップＳ４１８）。その後、着目する流体粒子２１について運動方程式を解き、その位置を時間発展させる（ステップＳ４１９）。

ステップＳ４１３において、着目する流体粒子２１に対して仮想粒子２５が対応付けられていると判定された場合、着目する流体粒子２１が仮想粒子除去条件を満たしているか否かを判定する（ステップＳ４１４）。

図８Ｂを参照して、仮想粒子除去条件について説明する。着目する流体粒子２１に対応する３個の仮想粒子２５が配置されている仮想的な平面２７から、着目する流体粒子２１までの距離ｒ_１が、斥力発生最長距離Ｌ_ｒｍ（図８Ｃ）以上である場合、着目する流体粒子２１は仮想粒子除去条件を満たしていると判定する。すなわち、仮想粒子除去条件が満たされる場合は、着目する流体粒子２１に、いずれの仮想粒子２５からも力が作用しない。なお、流体粒子２１と３個の仮想粒子２５との距離の最小値が斥力発生最長距離Ｌ_ｒｍ以上のとき、仮想粒子除去条件が満たされると判定するようにしてもよい。

ステップＳ４１４で仮想粒子除去条件が満たされていると判定された場合には、着目する流体粒子２１に対応付けられている仮想粒子２５（図８Ａ、図８Ｂ）を除去する（ステップＳ４１７）。仮想粒子２５を除去した状態で、着目する流体粒子２１に作用する力を計算する（ステップＳ４１８）。その後、着目する流体粒子２１について運動方程式を解き、その位置を時間発展させる（ステップＳ４１９）。

ステップＳ４１４で仮想粒子除去条件が満たされていないと判定された場合には、仮想粒子２５（図８Ａ、図８Ｂ）が配置されている条件の下で、着目する流体粒子２１に作用する力を計算する（ステップＳ４１６）。その後、着目する流体粒子２１について運動方程式を解き、その位置を時間発展させる（ステップＳ４１９）。

次に、上記実施例の優れた効果について説明する。
上記実施例では、壁面１１に近接した壁面近接粒子２１Ａ（図８Ａ、図８Ｂ）に対して複数の仮想粒子２５を配置し、仮想粒子２５が壁面近接粒子２１Ａに対して、壁面１１に平行な方向への移動を抑制する力を与えている。このため、解析モデルにおいて、壁面１１の近傍において流体の滑りなしの状態を再現することができる。

また、上記実施例では、壁面１１の形状に基づいて符号付距離関数を生成している。解析中に符号付距離関数を用いることにより、流体粒子２１から壁面１１までの距離や、流体粒子２１から壁面１１に下した垂線２６（図８Ｂ）の方向を容易に求めることができる。

次に図９を参照して、上記実施例の優れた効果を確認するために行ったシミュレーション結果について説明する。シミュレーションにおいて、円管１０の壁面１１（図３Ａ）を、半径２０ｍｍ、長さ８０ｍｍの円筒とした。長さ方向の両端面には周期境界条件を適用した。また、流体粒子２１に円管１０の軸方向の体積力を加え、平均流速が２ｍ／ｓになるまで解析を行った。高分子２０の挙動の解析には、クレマーグレストモデルを適用した。

図９は、べき乗則流体の流れをシミュレーションしたときのシミュレーション結果を示すグラフである。横軸は円管１０の中心軸からの半径方向の距離を単位「ｍｍ」で表し、縦軸は流体粒子２１の軸方向の速度を単位「ｍ／ｓ」で表す。図９に示したグラフにおいて、丸記号は実施例によるシミュレーション方法を用いて解析した結果を示し、ひし形記号は、比較例によるシミュレーション方法を用いて解析した結果を示し、実線は理論解を示す。

次に、比較例によるシミュレーション方法について簡単に説明する。比較例では、円管１０の壁面１１に沿って複数の壁面粒子を配置し、壁面粒子と流体粒子との間で相互作用ポテンシャルを定義する。１層目の複数の壁面粒子を、壁面１１に沿い、かつ流体粒子が壁面粒子の間をすり抜けることができる隙間を持たせて配置する。１層目の壁面粒子より深い位置に、１層目の壁面粒子の隙間を塞ぐように２層目の複数の壁面粒子を配置する。

流体粒子２１が１層目の複数の壁面粒子の隙間に入り込むことにより、壁面１１の近傍における流体の滑りなし状態が再現される。２層目の壁面粒子が１層目の壁面粒子の隙間を塞ぐことにより、流体粒子２１が壁面１１を貫通して円管１０の外側まで流出してしまうことが防止される。

実施例によるシミュレーション方法で求めた流速の分布は、実線で示した理論解によく一致していることがわかる。このシミュレーションから、実施例によるシミュレーション方法は、壁面近傍における流体の滑りなし状態をよく再現していることが確認された。

また、比較例による方法でも、壁面近傍における流体の滑りなし状態がよく再現されている。ところが、比較例によるシミュレーション方法では、壁面１１に沿う複数の壁面粒子を配置しなければならない。壁面粒子を配置するには、一般的に壁面１１を三角メッシュで区切り、節点に１層目の壁面粒子を配置する。生成される三角メッシュの品質によっては、２層目の壁面粒子で１層目の壁面粒子の隙間を埋めることが困難な場合が生じ得る。特に、壁面１１の幾何学的形状が複雑な場合に、三角メッシュの品質が低下しやすい。三角メッシュの品質が低下すると、１層目の壁面粒子の隙間を塞ぐように、２層目の壁面粒子を配置するためのアルゴリズムが複雑になる。

本実施例では、壁面１１に沿う壁面粒子を配置する必要がないため、２層目の壁面粒子を配置するための複雑なアルゴリズムが不要である。

次に、上記実施例の変形例について説明する。
上記実施例では、壁面１１（図３Ａ）を円筒形状としているが、壁面１１の形状は円筒形状に限らず、より複雑な形状の壁面についても、上記実施例を適用することが可能である。また、上記実施例では、壁面近接粒子２１Ａの近傍に３個の仮想粒子２５（図８Ａ）を配置したが、４個以上の仮想粒子２５を配置してもよい。

上記実施例では、壁面１１の形状を定義する符号付距離関数を利用して解析を行ったが、流体粒子２１から壁面１１までの距離や、流体粒子２１から壁面１１に下した垂線の方向を求めることができるその他の関数を用いて壁面１１の形状を定義してもよい。

上記実施例の効果を確認するためのシミュレーションでは、複数の流体粒子２１が高分子２０を構成しているべき乗則流体について解析を行ったが、上記実施例は、ニュートン流体の解析に適用することも可能である。

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

１０ 流体が流れる円管
１１ 壁面
２０ 高分子
２１ 流体粒子
２１Ａ 壁面近接粒子
２５ 仮想粒子
２６ 壁面近接粒子から壁面に下した垂線
２７ 仮想粒子が配置されている仮想的な平面
３０ 入力部
３１ 処理部
３２ 出力部
３３ 記憶部

Claims (9)

1. 壁面に接触する流体を複数の流体粒子で表した解析モデルにおいて、前記複数の流体粒子の挙動を解析するシミュレーション装置であって、
   前記壁面の形状を定義する情報、前記複数の流体粒子が前記壁面から受ける相互作用ポテンシャルを定義する情報、及び前記流体の物性値を含むシミュレーション条件が入力される入力部と、
   前記入力部に入力されたシミュレーション条件を取得し、取得された情報に基づいて前記複数の流体粒子について運動方程式を解いて前記複数の流体粒子の位置を時間発展させる処理部と
   を有し、
   前記処理部は、
   前記複数の流体粒子のうち前記壁面までの距離が近接判定閾値以下になった流体粒子を壁面近接粒子として検出し、
   前記壁面近接粒子と相互作用を及ぼし合う位置に複数の仮想粒子を発生させ、前記複数の仮想粒子の位置は固定し、前記壁面近接粒子と前記複数の仮想粒子との間に、前記壁面に対して平行な方向への前記壁面近接粒子の移動を妨げる相互作用ポテンシャルを作用させて、前記壁面近接粒子について運動方程式を解くシミュレーション装置。
2. 前記壁面近接粒子と前記複数の仮想粒子の各々との間に作用する相互作用ポテンシャルは、粒子間の距離が近づくにしたがって斥力が大きくなる形状を有する請求項１に記載のシミュレーション装置。
3. 前記処理部は、前記壁面近接粒子が前記複数の仮想粒子の各々から力を受けない距離まで遠ざかると、前記壁面近接粒子に対して発生した前記複数の仮想粒子を除去する請求項１または２に記載のシミュレーション装置。
4. 前記処理部は、前記複数の仮想粒子を発生させるときに、前記壁面近接粒子から前記壁面におろした垂線に対して直交する平面上であって、前記壁面近接粒子の位置を重心とする正三角形の頂点の位置にそれぞれ前記複数の仮想粒子を発生させる請求項１乃至３のいずれか１項に記載のシミュレーション装置。
5. 前記壁面近接粒子と前記複数の仮想粒子の各々との間に作用する相互作用ポテンシャルは、粒子間の距離が斥力発生最長距離以上のときは粒子に力を作用させず、粒子間の距離が前記斥力発生最長距離未満の時は粒子に斥力を作用させる請求項１乃至４のいずれか１項に記載のシミュレーション装置。
6. 前記処理部は、前記複数の仮想粒子を発生させるときに、前記壁面近接粒子からの距離が前記斥力発生最長距離の位置に前記複数の仮想粒子を発生させる請求項５に記載のシミュレーション装置。
7. 前記処理部は、前記複数の流体粒子が配置される空間を直交格子で分割し、前記壁面の形状を定義する情報に基づいて、各格子点に、前記壁面からの距離を対応付けた符号付距離関数を生成し、前記複数の流体粒子と前記壁面との距離を前記符号付距離関数を用いて求める請求項１乃至６のいずれか１項に記載のシミュレーション装置。
8. 壁面に接触する流体を複数の流体粒子で表した解析モデルにおいて、前記複数の流体粒子の挙動を解析するシミュレーション方法であって、
   前記壁面の形状を定義する情報、前記複数の流体粒子が前記壁面から受ける相互作用ポテンシャルを定義する情報、及び前記流体の物性値を含むシミュレーション条件を取得し、
   取得された情報に基づいて前記複数の流体粒子について運動方程式を解いて前記複数の流体粒子の挙動を解析し、
   解析中に、前記複数の流体粒子のうち前記壁面までの距離が近接判定閾値以下になった流体粒子を壁面近接粒子として検出し、
   検出された前記壁面近接粒子と相互作用を及ぼし合う位置に複数の仮想粒子を発生させ、前記複数の仮想粒子の位置は固定し、前記壁面近接粒子と前記複数の仮想粒子の間に、前記壁面近接粒子の前記壁面に対して平行な方向への移動を妨げる相互作用ポテンシャルを作用させて、前記壁面近接粒子について運動方程式を解くシミュレーション方法。
9. 壁面に接触する流体を複数の流体粒子で表した解析モデルにおいて、前記複数の流体粒子の挙動を解析するシミュレーションをコンピュータに実行させるプログラムであって、
   前記壁面の形状を定義する情報、前記複数の流体粒子が前記壁面から受ける相互作用ポテンシャルを定義する情報、及び前記流体の物性値を含むシミュレーション条件を取得する機能と、
   取得された情報に基づいて前記複数の流体粒子について運動方程式を解いて前記複数の流体粒子の挙動を解析する機能と、
   解析中に、前記複数の流体粒子のうち前記壁面までの距離が近接判定閾値以下になった流体粒子を壁面近接粒子として検出する機能と、
   検出された前記壁面近接粒子と相互作用を及ぼし合う位置に複数の仮想粒子を発生させ、前記複数の仮想粒子の位置は固定し、前記壁面近接粒子と前記複数の仮想粒子の間に、前記壁面近接粒子の前記壁面に対して平行な方向への移動を妨げる相互作用ポテンシャルを作用させて、前記壁面近接粒子について運動方程式を解く機能と
   をコンピュータに実現させるプログラム。
   

Images