特開 2023-175377 シミュレーション方法、シミュレーション装置、及びプログラム

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】公開特許公報(A)

(11)【公開番号】P2023175377

(43)【公開日】2023-12-12

(54)【発明の名称】シミュレーション方法、シミュレーション装置、及びプログラム

(51)【国際特許分類】

   G16Z 99/00 20190101AFI20231205BHJP

【ＦＩ】

G16Z99/00

【審査請求】未請求

【請求項の数】6

【出願形態】ＯＬ

(21)【出願番号】P 2022087797

(22)【出願日】2022-05-30

(71)【出願人】

【識別番号】000002107

【氏名又は名称】住友重機械工業株式会社

(74)【代理人】

【識別番号】100105887

【弁理士】

【氏名又は名称】来山 幹雄

(72)【発明者】

【氏名】廣瀬 良太

(72)【発明者】

【氏名】松宮 就章

【テーマコード（参考）】

5L049

【Ｆターム（参考）】

5L049DD02

(57)【要約】

【課題】壁面を再現するための壁粒子を配置することなく、壁面に接触する流体の流れを解析することができるシミュレーション方法を提供する。
【解決手段】壁面に接して流れる流体を複数の粒子で表す。複数の粒子と壁面との間の粒子壁面間相互作用、及び複数の粒子の間の粒子間相互作用を決定し、複数の粒子のそれぞれについて、複数の粒子の運動を支配する運動方程式を解くことにより、複数の粒子の位置及び速度を時間発展させる。運動方程式を解く際に、複数の粒子のうち壁面までの距離が、シミュレーション条件として設定された第１距離以下の粒子について、粒子間相互作用及び粒子壁面間相互作用による力の他に、壁面から受ける減衰力、及び壁面の温度に応じたランダム力を作用させて粒子の位置及び速度を時間発展させる。
【選択図】図４

【特許請求の範囲】

【請求項1】

壁面に接して流れる流体を複数の粒子で表し、
前記複数の粒子と前記壁面との間の粒子壁面間相互作用、及び前記複数の粒子の間の粒子間相互作用を決定し、
前記複数の粒子のそれぞれについて、前記複数の粒子の運動を支配する運動方程式を解くことにより、前記複数の粒子の位置及び速度を時間発展させるシミュレーション方法であって、
前記運動方程式を解く際に、
前記複数の粒子のうち前記壁面までの距離が、シミュレーション条件として設定された第１距離以下の粒子について、前記粒子間相互作用及び前記粒子壁面間相互作用による力の他に、前記壁面から受ける減衰力、及び前記壁面の温度に応じたランダム力を作用させて粒子の位置及び速度を時間発展させるシミュレーション方法。

【請求項2】

前記流体は一方向に流れ、
前記流体の流れの方向をｘ方向とするｘｙｚ直交座標系を定義したとき、前記複数の粒子は、ｘ方向、ｙ方向、及びｚ方向の少なくとも一方向にくりこまれたものであり、
ｘ方向、ｙ方向、ｚ方向へのくりこみ回数をそれぞれｎ_ｘ、ｎ_ｙ、ｎ_ｚと標記し、
くりこみの程度を表すくりこみ因子λ_ｘ、λ_ｙ、λ_ｚを、

【数1】

と標記し、
前記減衰力のくりこみ前の減衰係数をγ、くりこみ後の減衰係数をγ_Ｒと標記したとき、変換則

【数2】

を適用してくりこみ後の減衰係数γ_Ｒを計算し、
前記運動方程式を解く際に、くりこみ後の減衰係数γ_Ｒを用いる請求項１に記載のシミュレーション方法。

【請求項3】

壁面に沿って流れる流体の流れを解析するシミュレーション装置であって、
シミュレーション条件が入力される入力部と、
前記入力部に入力された前記シミュレーション条件に基づいて、前記流体の流れを解析する処理部と、
前記処理部による解析結果が出力される出力部と
を備え、
前記処理部は、前記入力部に入力された前記シミュレーション条件に基づいて前記流体を前記複数の粒子で表し、
前記複数の粒子のそれぞれについて、前記複数の粒子の運動を支配する運動方程式を解くことにより、前記複数の粒子の位置及び速度を時間発展させ、
前記運動方程式を解く際に、
前記複数の粒子のうち前記壁面までの距離が、前記シミュレーション条件で設定された第１距離以下の粒子について、前記シミュレーション条件で設定された粒子間相互作用及び粒子壁面間相互作用による力、前記壁面から受ける減衰力、及び前記壁面の温度に応じたランダム力を作用させて粒子の位置及び速度を時間発展させるシミュレーション装置。

【請求項4】

前記シミュレーション条件は、前記複数の粒子のくりこみを行うくりこみ条件を含み、
前記流体は一方向に流れ、
前記流体の流れの方向をｘ方向とするｘｙｚ直交座標系を定義したとき、前記くりこみ条件は、ｘ方向、ｙ方向、及びｚ方向へのくりこみ回数ｎ_ｘ、ｎ_ｙ、ｎ_ｚを指定する情報を含み、
くりこみの程度を表すくりこみ因子λ_ｘ、λ_ｙ、λ_ｚを、

【数3】

と標記し、
前記減衰力のくりこみ前の減衰係数をγ、くりこみ後の減衰係数をγ_Ｒと標記したとき、前記処理部は、変換則

【数4】

を適用してくりこみ後の減衰係数γ_Ｒを計算し、
前記運動方程式を解く際に、くりこみ後の減衰係数γ_Ｒを用いる請求項３に記載のシミュレーション装置。

【請求項5】

壁面に沿って流れる流体の流れを解析する手順をコンピュータに実行させるプログラムであって、
シミュレーション条件を取得する手順と、
取得された前記シミュレーション条件に基づいて、前記流体の流れを解析する手順と
を実行させ、
前記流体の流れを解析する手順は、
取得された前記シミュレーション条件に基づいて前記流体を前記複数の粒子で表す手順と、
前記複数の粒子のそれぞれについて、前記複数の粒子の運動を支配する運動方程式を解くことにより、前記複数の粒子の位置及び速度を時間発展させる手順と
を含み、
前記運動方程式を解く際に、
前記複数の粒子のうち前記壁面までの距離が、前記シミュレーション条件で設定される第１距離以下の粒子について、前記シミュレーション条件で設定された粒子間相互作用及び粒子壁面間相互作用による力、前記壁面から受ける減衰力、及び前記壁面の温度に応じたランダム力を作用させて粒子の位置及び速度を時間発展させるプログラム。

【請求項6】

前記シミュレーション条件は、前記複数の粒子のくりこみを行うくりこみ条件を含み、
前記流体は一方向に流れ、
前記流体の流れの方向をｘ方向とするｘｙｚ直交座標系を定義したとき、前記くりこみ条件は、ｘ方向、ｙ方向、及びｚ方向へのくりこみ回数ｎ_ｘ、ｎ_ｙ、ｎ_ｚを指定する情報を含み、
くりこみの程度を表すくりこみ因子λ_ｘ、λ_ｙ、λ_ｚを、

【数5】

と標記し、
前記減衰力のくりこみ前の減衰係数をγ、くりこみ後の減衰係数をγ_Ｒと標記したとき、
前記流体の流れを解析する手順は、
変換則

【数6】

を適用してくりこみ後の減衰係数γ_Ｒを計算する手順を含み、
前記運動方程式を解く際に、くりこみ後の減衰係数γ_Ｒを用いる請求項５に記載のプログラム。

【発明の詳細な説明】

【技術分野】

【0001】

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

【背景技術】

【0002】

壁面に接する流体の流れを、分子動力学法を用いて解析する従来の方法では、壁面を複数の壁粒子で表して、壁粒子と流体粒子との間の相互作用を考慮する（特許文献１参照）。壁粒子の温度制御を行うことにより、壁面近傍の流体粒子の温度が再現される。

【先行技術文献】

【特許文献】

【0003】

【特許文献1】特開２００７―２１９８３１号公報

【発明の概要】

【発明が解決しようとする課題】

【0004】

従来のシミュレーション方法では、流体粒子の他に壁粒子を配置する必要があるため、計算対象の粒子数が多くなる。また、壁面形状の空間分解能が、壁粒子の粒子サイズによって制限される。

【0005】

本発明の目的は、壁面を再現するための壁粒子を配置することなく、壁面に接触する流体の流れを解析することができるシミュレーション方法、シミュレーション装置、及びプログラムを提供することである。

【課題を解決するための手段】

【0006】

本発明の一観点によると、
壁面に接して流れる流体を複数の粒子で表し、
前記複数の粒子と前記壁面との間の粒子壁面間相互作用、及び前記複数の粒子の間の粒子間相互作用を決定し、
前記複数の粒子のそれぞれについて、前記複数の粒子の運動を支配する運動方程式を解くことにより、前記複数の粒子の位置及び速度を時間発展させるシミュレーション方法であって、
前記運動方程式を解く際に、
前記複数の粒子のうち前記壁面までの距離が、シミュレーション条件として設定された第１距離以下の粒子について、前記粒子間相互作用及び前記粒子壁面間相互作用による力の他に、前記壁面から受ける減衰力、及び前記壁面の温度に応じたランダム力を作用させて粒子の位置及び速度を時間発展させるシミュレーション方法が提供される。

【0007】

本発明の他の観点によると、
壁面に沿って流れる流体の流れを解析するシミュレーション装置であって、
シミュレーション条件が入力される入力部と、
前記入力部に入力された前記シミュレーション条件に基づいて、前記流体の流れを解析する処理部と、
前記処理部による解析結果が出力される出力部と
を備え、
前記処理部は、前記入力部に入力された前記シミュレーション条件に基づいて前記流体を前記複数の粒子で表し、
前記複数の粒子のそれぞれについて、前記複数の粒子の運動を支配する運動方程式を解くことにより、前記複数の粒子の位置及び速度を時間発展させ、
前記運動方程式を解く際に、
前記複数の粒子のうち前記壁面までの距離が、前記シミュレーション条件で設定された第１距離以下の粒子について、前記シミュレーション条件で設定された粒子間相互作用及び粒子壁面間相互作用による力、前記壁面から受ける減衰力、及び前記壁面の温度に応じたランダム力を作用させて粒子の位置及び速度を時間発展させるシミュレーション装置が提供される。

【0008】

本発明のさらに他の観点によると、
壁面に沿って流れる流体の流れを解析する手順をコンピュータに実行させるプログラムであって、
シミュレーション条件を取得する手順と、
取得された前記シミュレーション条件に基づいて、前記流体の流れを解析する手順と
を実行させ、
前記流体の流れを解析する手順は、
取得された前記シミュレーション条件に基づいて前記流体を前記複数の粒子で表す手順と、
前記複数の粒子のそれぞれについて、前記複数の粒子の運動を支配する運動方程式を解くことにより、前記複数の粒子の位置及び速度を時間発展させる手順と
を含み、
前記運動方程式を解く際に、
前記複数の粒子のうち前記壁面までの距離が、前記シミュレーション条件で設定される第１距離以下の粒子について、前記シミュレーション条件で設定された粒子間相互作用及び粒子壁面間相互作用による力、前記壁面から受ける減衰力、及び前記壁面の温度に応じたランダム力を作用させて粒子の位置及び速度を時間発展させるプログラムが提供される。

【発明の効果】

【0009】

壁面を再現する壁面粒子を配置しないため、計算対象の粒子数の増大が抑制され、計算負荷を低減させることができる。流体の粒子に、壁面から受ける減衰力、及び壁面の温度に応じたランダム力を作用させることにより、壁面におけるノンスリップ条件を再現し、流体の粒子の温度を壁面の温度に維持することができる。

【図面の簡単な説明】

【0010】

【図1】図１Ａ及び図１Ｂは、一実施例によるシミュレーション方法で解析対象となる解析モデルの一例を模式的に示す断面図である。

【図2】図２Ａは、壁面からの距離が第１距離より遠い位置の粒子に作用する力を示す模式図であり、図２Ｂは、壁面からの距離が第１距離以下の位置にある粒子に作用する力を示す模式図である。

【図3】図３は、本実施例によるシミュレーション装置のブロック図である。

【図4】図４は、本実施例によるシミュレーション方法の手順を示すフローチャートである。

【図5】図５Ａは、解析対象の流体を複数の粒子で表した解析モデルの模式図であり、図５Ｂは、図５Ａに示した粒子系に対して、等方的なくりこみを行った粒子系を示す模式図であり、図５Ｃは、図５Ｂに示した粒子系に対して、ｙ方向にはくりこみを行わず、ｘ方向及びｚｚ方向にさらにくりこみを行った粒子系を示す模式図である。

【図6】図６は、くりこみを行わない解析モデルを用いた解析結果から計算したｘ方向の流速のｙ方向の分布を示すグラフである。

【図7】図７は、ｘ、ｙ、ｚ方向に１回のくりこみを行った解析モデルを用いた解析結果から計算したｘ方向の流速のｙ方向の分布を示すグラフである。

【図8】図８は、ｘ方向に１回のくりこみを行った解析モデルを用いた解析結果から計算したｘ方向の流速のｙ方向の分布を示すグラフである。

【図9】図９は、ｙ方向に１回のくりこみを行った解析モデルを用いた解析結果から計算したｘ方向の流速のｙ方向の分布を示すグラフである。

【図10】図１０は、ｚ方向に１回のくりこみを行った解析モデルを用いた解析結果から計算したｘ方向の流速のｙ方向の分布を示すグラフである。

【図11】図１１は、ｘ、ｙ、ｚ方向に２回のくりこみを行った解析モデルを用いた解析結果から計算したｘ方向の流速のｙ方向の分布を示すグラフである。

【図12】図１２は、ｘ方向に２回のくりこみを行った解析モデルを用いた解析結果から計算したｘ方向の流速のｙ方向の分布を示すグラフである。

【図13】図１３は、ｙ方向に２回のくりこみを行った解析モデルを用いた解析結果から計算したｘ方向の流速のｙ方向の分布を示すグラフである。

【図14】図１４は、ｚ方向に２回のくりこみを行った解析モデルを用いた解析結果から計算したｘ方向の流速のｙ方向の分布を示すグラフである。

【図15】図１５は、ｘ、ｙ、ｚ方向に３回のくりこみを行った解析モデルを用いた解析結果から計算したｘ方向の流速のｙ方向の分布を示すグラフである。

【図16】図１６は、ｘ方向に３回のくりこみを行った解析モデルを用いた解析結果から計算したｘ方向の流速のｙ方向の分布を示すグラフである。

【図17】図１７は、ｙ方向に３回のくりこみを行った解析モデルを用いた解析結果から計算したｘ方向の流速のｙ方向の分布を示すグラフである。

【図18】図１８は、ｚ方向に３回のくりこみを行った解析モデルを用いた解析結果から計算したｘ方向の流速のｙ方向の分布を示すグラフである。

【図19】図１９は、ｘ、ｙ、ｚ方向に３回のくりこみを行った解析モデルを用いた解析結果から計算したｘ方向の流速のｙ方向の分布を示すグラフである。

【図20】図２０は、ｘ方向に３回のくりこみを行った解析モデルを用い、比較例による方法で解析を行った結果から計算したｘ方向の流速のｙ方向の分布を示すグラフである。

【図21】図２１は、ｙ方向に３回のくりこみを行った解析モデルを用い、比較例による方法で解析を行った結果から計算したｘ方向の流速のｙ方向の分布を示すグラフである。

【図22】図２２は、ｚ方向に３回のくりこみを行った解析モデルを用い、比較例による方法で解析を行った結果から計算したｘ方向の流速のｙ方向の分布を示すグラフである。

【図23】図２３Ａ及び図２３Ｂは、他の実施例によるシミュレーション方法の解析対象となる解析モデルの断面図である。

【発明を実施するための形態】

【0011】

図１Ａから図４までの図面を参照して、一実施例によるシミュレーション方法、シミュレーション装置、及びプログラムについて説明する。

【0012】

図１Ａ及び図１Ｂは、一実施例によるシミュレーション方法の解析対象となる解析モデルの一例を模式的に示す断面図である。図１Ａは、図１Ｂの一点鎖線１Ａ－１Ａにおける断面図であり、図１Ｂは、図１Ａの一点鎖線１Ｂ－１Ｂにおける断面図である。壁面４０で画定される解析空間３０内を流体が流れる。壁面４０は、例えば相互に平行に配置された一対の面で構成される。壁面４０に平行な面をｘｚ面とし、流体の流れの方向をｘ方向とするｘｙｚ直交座標系を定義する。

【0013】

壁面４０の間の解析空間３０のｘ、ｙ、ｚ方向の寸法を、それぞれＬｘ、Ｌｙ、Ｌｚと標記する。ｚ方向に垂直な境界、及びｘ方向に垂直な境界には、周期境界条件を適用する。解析空間３０内の流体を複数の粒子３１で表す。実施例においては、分子動力学法またはくりこみ群分子動力学法を用いて、複数の粒子３１の挙動を解析する。具体的には、複数の粒子３１の運動を支配する運動方程式を数値的に解いて、粒子３１の位置及び速度を時間発展させる。

【0014】

粒子３１の間に作用する相互作用ポテンシャルＵとして、例えば以下のレナードジョーンズポテンシャルを用いる。

【数1】

ここで、ｒは粒子間の距離であり、ε及びσは、それぞれエネルギ及び長さの次元を持つフィッティングパラメータである。σは、粒子の衝突直径といわれる場合がある。

【0015】

図２Ａは、壁面４０からの距離Ｌ_ｉｗが第１距離Ｌ_１より遠い位置の粒子３１に作用する力を示す模式図である。ここで、添え字のｉは、複数の粒子３１に通し番号を付したときのｉ番目の粒子３１を意味する。

【0016】

壁面４０からの距離Ｌ_ｉｗが第１距離Ｌ_１より遠い位置の粒子３１の運動を支配する運動方程式として、以下の式を適用することができる。

【数2】

ここで、ｍは粒子３１の質量であり、ｖ_ｉは、ｉ番目の粒子３１の速度ベクトルである。Ｆ_ｉｊは、ｉ番目の粒子３１がｊ番目の粒子３１から受ける力であり、式（１）に示したレナードジョーンズポテンシャルから計算することができる。Ｆ_ｅｘｔは、粒子３１が受ける重力等の外力である。式（２）の右辺のΣ記号は、ｉ番目の粒子３１の周囲の粒子３１についての合計を意味する。

【0017】

壁面４０を複数の壁粒子で表して、粒子３１と壁粒子との間の相互作用を考慮することにより、壁面で流速がゼロになるノンスリップ条件が再現される。本実施例では、壁面を壁粒子で表さず、ノンスリップ条件を再現する。

【0018】

図２Ｂは、壁面４０からの距離Ｌ_ｉｗが第１距離Ｌ_１以下の位置にある粒子３１（以下、単に、壁面近傍の粒子３１という場合がある。）に作用する力を示す模式図である。ノンスリップ条件を再現するために、壁面近傍の粒子３１を支配する運動方程式に減衰項を追加すると、時間の経過とともに壁面近傍の粒子３１の温度が低下してしまう。本実施例では、ランジュバン（Langevin）法を用いて壁面近傍の粒子３１の温度制御を行う。

【0019】

ランジュバン法による温度制御を適用した場合に、壁面近傍の粒子３１の運動を支配する運動方程式は、以下の式で表される。

【数3】

ここで、Ｆ_ｉｗは、ｉ番目の粒子３１が壁面４０から受ける法線方向の力である。ｖ_ｉ、ｖ_ｗは、それぞれｉ番目の粒子３１の速度ベクトル、及びｉ番目の粒子３１が接近している壁面４０の速度ベクトルである。式（３）の右辺第４項は減衰力（粘性抵抗力）に相当する。減衰力の大きさは、壁面４０に対する粒子３１の相対速度に比例し、減衰力の向きは、相対速度と反対向きである。γを減衰係数ということとする。この減衰力により、壁面４０におけるノンスリップ条件を再現することができる。

【0020】

Ｒ_ｉは、壁面４０の近傍の粒子３１の温度を、壁面温度に維持するために付与するランダム力である。ランダム力Ｒ_ｉの大きさは、平均がゼロ、標準偏差σ_ｓが以下の式で表される正規分布に従う。

【数4】

ここで、ｋ_Ｂはボルツマン定数であり、Ｔ_ｗは壁面４０の設定温度であり、Δｔ_ｓは運動方程式を数値的に解くときの時間刻み幅である。ランダム力Ｒ_ｉの向きは、ランダムである。

【0021】

ランジュバン法による温度制御は、減衰係数γによって制御の強さが変化する。減衰係数γを小さく設定すると温度制御が弱くなり、粒子３１の温度の追従性が悪くなる。逆に減衰係数γを大きくすると、時間刻み幅Δｔ_ｓとの兼ね合いで、設定温度Ｔ_ｗからずれが生じたり、計算が破綻したりする場合がある。

【0022】

減衰係数γは、例えば以下の式で表される。

【数5】

ここで、ζは減衰比であり、ｋはバネ定数である。減衰比ζとして、例えば０．７０７を用いることができる。バネ定数ｋは、例えば、レナードジョーンズポテンシャル（式（１））を鞍点近傍でテイラー展開したときの２次の項の係数から求めることができる。この場合、バネ定数ｋは、以下の式で計算することができる。

【数6】

減衰係数γとして、例えば５．２３×１０－１３ｋｇ／ｓを用いることができる。

【0023】

次に、ｉ番目の粒子３１が壁面４０から受ける法線方向の力Ｆ_ｉｗについて説明する。粒子３１から壁面４０までの距離Ｌ_ｉｗが第１距離Ｌ_１以下になると、壁面４０に関して面対称の位置に仮想的な粒子３１Ｖを配置する。第１距離Ｌ_１として、例えば式（１）の衝突直径σの１／２を採用する。なお、図２Ｂにおいて、粒子３１を表す円形の直径は、粒子３１の衝突直径σを表しているわけではない。

【0024】

壁面近傍のｉ番目の粒子３１に、仮想的な粒子３１Ｖから、式（１）で定義される粒子間相互作用ポテンシャルに基づく力Ｆ_ｉｗを作用させる。すなわち、粒子３１が衝突直径σの１／２以下となる距離まで壁面４０に近づくと、壁面４０から法線方向の斥力を受ける。

【0025】

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

【0026】

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

【0027】

次に、図４を参照して本実施例によるシミュレーション方法について説明する。図４は、本実施例によるシミュレーション方法の手順を示すフローチャートである。フローチャートに示された各手順は、処理部５１が記憶部５３に記憶されているプログラムを実行することにより行われる。

【0028】

まず、処理部５１が入力部５０に入力されたシミュレーション条件を取得する（ステップＳ１）。シミュレーション条件には、解析対象の流体を定義する情報、壁面４０（図１Ａ、図１Ｂ）の形状を定義する情報、流体を複数の粒子３１（図１Ａ、図１Ｂ）で表すための情報、複数の粒子３１と壁面４０との間の粒子壁面間相互作用を定義する情報、複数の粒子３１の間の粒子間相互作用を定義する情報、温度条件を定義する情報、壁面近傍の粒子３１が壁面４０から受ける減衰力を定義する情報、粒子に作用する外力を定義する情報、時間刻み幅、解析終了条件等が含まれる。

【0029】

処理部５１は、入力されたシミュレーション条件に基づいて解析空間３０（図１Ａ、図１Ｂ）内に複数の粒子３１を配置する。その後、ステップＳ４、ステップＳ５、ステップＳ６の手順を、すべての粒子３１について繰り返す（ステップＳ３）。以下、ステップＳ３の繰返し処理を、ｉ番目の粒子３１に着目して説明する。

【0030】

まず、ｉ番目の粒子３１から壁面４０までの距離Ｌ_ｉｗが第１距離Ｌ_１以下か否かを判定する（ステップＳ４）。ｉ番目の粒子３１から壁面４０までの距離Ｌ_ｉｗが第１距離Ｌ_１より長い場合、図２Ａに示した粒子間相互作用による力Ｆ_ｉｊ及び外力Ｆ_ｅｘｔを考慮して、粒子３１の運動を支配する運動方程式（式（２））を数値的に解く（ステップＳ５）。ｉ番目の粒子３１から壁面４０までの距離Ｌ_ｉｗが第１距離Ｌ_１以下の場合、図２Ｂに示した粒子間相互作用による力Ｆ_ｉｊ、粒子壁面間相互作用による力Ｆ_ｉｗ、粒子３１と壁面との相対速度に応じた減衰力－γ（ｖ_ｉ－ｖ_ｗ）、壁面４０の設定温度に応じたランダム力Ｒ_ｉ、及び外力Ｆ_ｅｘｔを考慮して、粒子３１の運動を支配する運動方程式（式（３））を数値的に解く（ステップＳ６）。

【0031】

すべての粒子３１について運動方程式を解いた後、粒子３１のそれぞれの位置及び速度を時間発展させる（ステップＳ７）。ステップＳ３の繰返し処理、及びステップＳ７を、解析が終了するまで繰り返す（ステップＳ８）。解析終了の条件は、例えば、ステップＳ１で取得するシミュレーション条件で与えられている。解析が終了すると、解析結果を出力部５２に出力する（ステップＳ９）。

【0032】

次に、図１Ａ～図４に示した実施例の優れた効果について説明する。上記実施例では、壁面４０（図１Ａ、図１Ｂ）を複数の壁粒子で表すことなく解析を行う。このため、解析対象の粒子数が少なくなる。その結果、計算負荷を軽減することができる。

【0033】

また、壁面近傍の粒子３１に壁面４０からの減衰力－γ（ｖ_ｉ－ｖ_ｗ）を作用させる（ステップＳ６、図２Ｂ）ことにより、壁面４０の表面におけるノンスリップ条件を再現することができる。さらに、壁面近傍の粒子３１に壁面４０の温度に応じたランダム力Ｒ_ｉを作用させる（ステップＳ６、図２Ｂ）ことにより、壁面近傍の粒子３１の温度を設定温度に維持することができる。

【0034】

上記実施例では、平行に配置された２枚の壁面内を流れる流体の流れについて解析を行っているが、壁面が他の形状を有する場合でも、上記実施例によるシミュレーション方法及びシミュレーション装置を適用して、流体の流れを解析することができる。

【0035】

次に、図５Ａ～図５Ｃを参照して他の実施例によるシミュレーション方法、シミュレーション装置、及びプログラムについて説明する。図１Ａ～図４を参照して説明した実施例では、粒子３１のくりこみを行っていない。以下に説明する実施例では、粒子３１のくりこみを行って粒子数を削減し、くりこみ後の粒子系に対して、式（３）及び式（４）に示した運動方程式を適用する。

【0036】

まず、本実施例で適用されるくりこみの手法について説明する。
図５Ａは、解析対象の流体を複数の粒子３１で表した解析モデルの模式図である。平行に配置された一対の壁面４０で挟まれた空間に、解析対象の流体が収容されている。流体は、複数の粒子３１で表される。粒子３１のそれぞれの形状は、粒子３１によって発生する相互作用ポテンシャルの等ポテンシャル面に相当すると考えることができる。一般的に、相互作用ポテンシャルの等ポテンシャル面の形状は球面形状である。図５Ａでは、ある大きさの等ポテンシャル面で粒子３１の各々を表している。壁面４０が隔たる方向をｙ方向とするｘｙｚ直交座標系を定義する。流体が収容される空間のｙ方向の寸法が、ｘ方向及びｚ方向の寸法に比べて十分小さい。

【0037】

図５Ｂは、図５Ａに示した粒子系に対して、等方的なくりこみを行った粒子系を示す模式図である。くりこみにより粒子数が減少する。相互作用ポテンシャルは、ｘ方向、ｙ方向、及びｚ方向に等方的に引き延ばされる。このため、図５Ｂにおいて、くりこみ後の粒子３２の各々を、図５Ａに示した粒子３１より大きな球面で表している。

【0038】

図５Ｃは、図５Ｂに示した粒子系に対して、ｙ方向にはくりこみを行わず、ｘ方向及びｚ方向にさらにくりこみを行った粒子系を示す模式図である。くりこみにより、粒子数がさらに減少する。また、相互作用ポテンシャルは、ｘ方向及びｚ方向に引き延ばされ、ｙ方向には引き延ばされない。このため、図５Ｃにおいて、くりこみ後の粒子３３の各々を、ｙ方向を短軸方向とする扁平な回転楕円体で表している。

【0039】

［等方的なくりこみ］
次に、粒子系に対して等方的なくりこみを行う場合のくりこみ変換則、及び粒子系のエネルギについて説明する。くりこみ変換則は、以下の式で表される。

【数7】

ここで、Ｎは粒子数であり、ｍは粒子の質量であり、ｒは粒子の位置を示す位置ベクトルであり、Ｖは流体の流れ場の体積であり、Ｔは粒子系の温度である。下付きの添え字Ｒが付されたパラメータは、くりこみ後のパラメータであることを表している。λは、くりこみの程度を表すパラメータ（くりこみ因子）であり、λは１より大きい実数である。例えば、ｎを１以上の整数として、くりこみ因子λを以下の式で表した時、ｎはくりこみ回数と呼ばれる。

【数8】

【0040】

粒子間距離をｒとして粒子間の相互作用ポテンシャルｕ（ｒ）が与えられる。距離ｒが無限大に近づくとき、ｕ（ｒ）が十分速くゼロに近づく場合（例えば、レナードジョーンズポテンシャルの場合）、相互作用ポテンシャルｕ（ｒ）のくりこみ変換則は、以下の式で表される。

【数9】

【0041】

上記の式（７）、式（９）のくりこみ変換則を用いてくりこみ処理を行った場合、くりこみの前後で粒子系のエネルギ、圧力等の巨視的な物理量は不変である。以下、くりこみの前後でエネルギが不変であることについて証明する。

【0042】

粒子系の分配関数Ｚ_Ｎは以下の式で表される。

【数10】

ここで、ｈはプランク定数、ｋ_Ｂはボルツマン定数、ｒは粒子の位置ベクトル、ｐは粒子の運動量、ｒ_ｉｊはｉ番目の粒子からｊ番目の粒子までの距離である。式（１０）のうち３番目の式の右辺の第１項のシグマは、Ｎ個のすべての粒子について合計することを意味しており、第２項のシグマは、すべての粒子対について合計することを意味している。

【0043】

式（１０）の分配関数Ｚ_Ｎの運動項の部分Ｚ_Ｎ：ｋは、以下の式で表される。

【数11】

【0044】

従って、運動エネルギＥ_ｋは、以下の式で表される。

【数12】

くりこみ因子λのくりこみを行うと、式（１２）の右辺のＮ及びβが共に１／λ^３になるため、エネルギＥ_Ｋは、くりこみの前後で不変である。

【0045】

式（１０）の分配関数Ｚ_Ｎの相互作用項の部分Ｚ_{Ｎ：ｉｎｔ}は、以下の式で表される。

【数13】

距離ｒ_ｉｊが無限大に近づくと、ｕ（ｒ_ｉｊ）は十分速くゼロに近づくため、あるカットオフ距離ｒ_ｃを用いて、ｒ＞ｒ_ｃのとき相互作用ポテンシャルｕ（ｒ）を以下のように近似することができる。

【数14】

【0046】

カットオフ距離ｒ_ｃが、流体を収容する空間のｘ、ｙ、ｚ方向の長さＬ_ｘ、Ｌ_ｙ、Ｌ_ｚより十分小さく、かつ粒子数密度Ｎ／Ｖが十分小さい場合を考える。このとき、各粒子からの距離がカットオフ距離ｒ_ｃ以下の範囲に他の粒子が存在する確率は極めて小さい。特に、各粒子からの距離がカットオフ距離ｒ_ｃ以下の範囲に２個以上の他の粒子が存在する確率はゼロとみなしてよい。したがって、分配関数Ｚ_{Ｎ：ｉｎｔ}の多重積分（式（１３））は、以下のように近似することができる。

【0047】

【数15】

【0048】

従って、相互作用エネルギＥ_ｉｎｔは、以下の式で記事することができる。

【数16】

【0049】

粒子系の全エネルギＥは、以下の式で定義される。

【数17】

カットオフ距離ｒ_ｃがＬ_ｘ、Ｌ_ｙ、Ｌ_ｚより十分小さく、粒子間距離ｒがカットオフ距離ｒ_ｃより大きいときｕ（ｒ）＝０と近似できるため、式（１６）の体積積分は、以下のように近似することができる。

【数18】

【0050】

くりこみ変換後の相互作用エネルギＥ_{ｉｎｔ，Ｒ}は、以下の式で近似される。

【数19】

【0051】

なお、式（１８）の近似が成立するためには、くりこみ変換後のカットオフ距離ｒ_ｃＲがＬ_ｘ、Ｌ_ｙ、Ｌ_ｚより十分小さいという条件を満たす必要がある。式（９）から、くりこみ変換後のカットオフ距離ｒ_ｃＲは、くりこみ変換前のカットオフ距離ｒ_ｃのλ倍である。くりこみ因子λを大きくすると、カットオフ距離ｒ_ｃＲが長くなる。このため、くりこみ因子λの上限は、くりこみ変換後のカットオフ距離ｒ_ｃＲがＬ_ｘ、Ｌ_ｙ、Ｌ_ｚより十分小さいという条件によって制約を受ける。

【0052】

式（１９）においてｒをλｒに変数変換すると、式（１９）は、式（１６）の右辺と同一になる。したがって、Ｅ_{ｉｎｔ，Ｒ}＝Ｅ_ｉｎｔが成立し、相互作用エネルギＥ_ｉｎｔも、くりこみの前後で不変である。

【0053】

このように、式（７）に示したくりこみ変換則を用いて等方的なくりこみを行うと、くりこみの前後で系の運動エネルギ及び相互作用エネルギは不変である。したがって、式（１７）に示した系全体のエネルギも、くりこみの前後で不変である。

【0054】

［二方向へのくりこみ］
次に、流れ場のｘ方向の寸法Ｌ_ｘ及びｚ方向の寸法Ｌ_ｚと比べてカットオフ距離ｒ_ｃが十分短いが、ｙ方向の寸法Ｌ_ｙと比べて十分短いとはいえない場合について説明する。例えば、流れ場が、ｙ方向を厚さ方向とする薄板状の形状である場合に、この条件が満たされる。

【0055】

流れ場の形状がこのような条件を満たす場合、ｙ方向にはくりこみを行わず、ｘ方向及びｚ方向の二方向にのみくりこみを行う。くりこみ変換則は以下の式で表される。

【数20】

【0056】

相互作用ポテンシャルは、以下の変換則に従う。

【数21】

【0057】

式（２０）のくりこみ変換を行う場合でも、式（１２）で表される運動エネルギはくりこみの前後で不変である。

【0058】

式（１６）の体積積分は、以下のように近似することができる。

【数22】

【0059】

くりこみ変換後の相互作用エネルギＥ_{ｉｎｔ，Ｒ}は、以下の式で近似される。

【数23】

ここで、ｒチルダは、式（２１）のｒチルダと同一である。式（２３）においてｘをλｘに変数変換し、ｚをλｚに変数変換すると、式（２３）の右辺は、式（１６）の右辺と同一になる。したがって、Ｅ_{ｉｎｔ，Ｒ}＝Ｅ_ｉｎｔが成立し、相互作用エネルギＥ_ｉｎｔも、くりこみの前後で不変である。

【0060】

このように、式（２０）に示した二方向についてのくりこみ変換則に基づいて異方的なくりこみを行うことにより、くりこみの前後で粒子系のエネルギを不変にすることができる。

【0061】

［一方向へのくりこみ］
次に、流れ場のｘ方向の寸法Ｌ_ｘと比べてカットオフ距離ｒ_ｃが十分短いが、ｙ方向の寸法Ｌ_ｙ及びｚ方向の寸法Ｌ_ｚと比べてカットオフ距離ｒ_ｃが十分短いとはいえない場合について説明する。例えば、流れ場が、ｘ方向を長さ方向とする細長い円柱形状である場合に、この条件が満たされる。

【0062】

流れ場の形状がこのような条件を満たす場合、ｙ方向及びｚ方向にはくりこみを行わず、ｘ方向にのみくりこみを行う。くりこみ変換則は以下の式で表される。

【数24】

【0063】

相互作用ポテンシャルは、以下の変換則に従う。

【数25】

【0064】

式（２４）のくりこみ変換を行う場合でも、式（１２）で表される運動エネルギはくりこみの前後で不変である。

【0065】

式（１６）の体積積分は、以下のように近似することができる。

【数26】

【0066】

くりこみ変換後の相互作用エネルギＥ_{ｉｎｔ，Ｒ}は、以下の式で近似される。

【数27】

ここで、ｒチルダは、式（２５）のｒチルダと同一である。式（２７）においてｘをλｘに変数変換すると、式（２７）の右辺は、式（１６）の右辺と同一になる。したがって、Ｅ_{ｉｎｔ，Ｒ}＝Ｅ_ｉｎｔが成立し、相互作用エネルギＥ_ｉｎｔも、くりこみの前後で不変である。

【0067】

このように、式（２４）に示した一方向についてのくりこみ変換則に基づいて異方的なくりこみを行うことにより、くりこみの前後で粒子系のエネルギを不変にすることができる。

【0068】

［異方的なくりこみの一般化］
二方向にくりこみを行い、残り一方向にはくりこみを行わない場合のくりこみ変換則を式（２０）に示し、二方向にはくりこみを行わず、残りの一方向にのみくりこみを行う場合のくりこみ変換則を式（２４）に示している。次に、三方向のそれぞれの方向に対してくりこみの程度を決定し、異方的なくりこみを行う場合について説明する。

【0069】

ｘ方向、ｙ方向、及びｚ方向へのくりこみ因子λ_ｘ、λ_ｙ、及びλ_ｚと標記する。このときのくりこみ変換則は、以下の式で表される。

【数28】

【0070】

相互作用ポテンシャルは、以下の変換則に従う。

【数29】

【0071】

くりこみ変換後の相互作用ポテンシャルｕ_Ｒ（ｒ）のｘ、ｙ、ｚ方向のカットオフ距離ｒ_ｃｘＲ、ｒ_ｃｙＲ、ｒ_ｃｚＲは、以下の式で表される。

【数30】

【0072】

このとき、式（１８）、式（２２）、式（２６）の近似と同様の近似が成立するようにするために、以下の条件を満たすことが好ましい。

【数31】

【0073】

次に、くりこみ後の粒子系に対して式（３）及び式（４）に示した運動方程式を適用するためには、減衰係数γの変換則を定義しなければならない。くりこみ後の減衰係数をγ_Ｒと標記する。等方的なくりこみを行う場合、以下の変換則を適用すると、くりこみ前の粒子系とくりこみ後の粒子系との間で相似な結果が得られる。

【数32】

【0074】

異方的にくりこみを行う場合、式（３２）から、以下の変換則を適用することが考えられる。

【数33】

【0075】

ところが、くりこみを行った粒子系について式（３３）の変換則を用いて得られた減衰係数γ_Ｒを用いてシミュレーションを行ったところ、くりこみの前後の粒子系の間で相似な結果が得られなかった。本願の発明者らによる種々の解析により、流体の流れの方向がｘ方向である場合に以下の変換則を適用すると、くりこみの前後で相似な結果が得られることがわかった。

【数34】

【0076】

次に、本実施例の優れた効果について説明する。
本実施例では、くりこみを行うことにより、粒子数を削減し、計算負荷を低減させることができる。さらに、解析空間の形状に応じて異方的なくりこみを行うことにより、さらに粒子数を削減し、計算負荷を低減させることができる。

【0077】

次に、図６～図２２を参照して、実際にシミュレーションを行った結果について説明する。図１Ａ、図１Ｂに示したポアズイユ流れの解析モデルについて、くりこみを行わない方法、等方的なくりこみを行う方法、異方的なくりこみを行う方法で解析を行った。ｘ方向とｚ方向とに、周期境界条件を適用した。ｘ方向の寸法Ｌｘ及びｚ方向の寸法Ｌｚをともに７．２０ｎｍとし、ｙ方向の寸法Ｌｙを６．２９ｎｍとした。

【0078】

くりこみ後も、無次元化された長さを同一に保つように解析モデルの作成を行った。解析対象の流体として水（Ｈ_２Ｏ）を想定し、式（１）のフィッティングパラメータεを４０４．５Ｋとし、σを０．２６４ｎｍとした。また、式（２）、式（３）の外力Ｆ_ｅｘｔとして、ｘ方向の外力を作用させた。外力の大きさは、流速の最大値が約１００ｍ／ｓになるように設定した。壁面４０の温度は、粒子系の温度の初期値と等しくした。

【0079】

流れが定常状態になるまで解析を行い、定常状態の粒子３１の移動速度からｘ方向の流速を計算した。図６～図２２は、解析結果から計算したｘ方向の流速の、ｙ方向に関する分布を示すグラフである。横軸はｙ方向の位置、すなわち一方の壁面４０からの距離を単位［ｎｍ］で表し、縦軸は流速のｘ方向成分を単位［ｍ／ｓ］で表す。各グラフの実線は理論値を示しており、丸記号は解析結果から求めた流速を示している。各グラフに付したｎ_ｘ、ｎ_ｙ、ｎ_ｚは、それぞれｘ方向、ｙ方向、ｚ方向のくりこみ回数を表す。すなわち、くりこみ回数は以下の式で定義される。

【数35】

【0080】

くりこみ回数ｎ_ｘ、ｎ_ｙ、ｎ_ｚ、またはくりこみ因子λ_ｘ、λ_ｙ、λ_ｚ等で指定されるくりこみ条件は、例えば図４のステップＳ１で取得するシミュレーション条件で与えられる。

【0081】

図６は、くりこみを行っていない場合の解析結果を示す。減衰係数γを、５．２３×１０^－１３ｋｇ／ｓとした。理論値に対する解析結果の平均誤差は約３．８％であり、解析結果は理論値とよく一致していることがわかる。図７は、等方的に１回ずつくりこみを行った場合の解析結果を示す。図８、図９、図１０は、それぞれｘ方向、ｙ方向、ｚ方向に１回のくりこみを行い、他の方向にはくりこみを行っていない場合の解析結果を示す。なお、図７～図１０では、式（３４）のくりこみ変換則を用いている。図７、図８、図９、図１０に示した解析結果の理論値に対する平均誤差は、それぞれ約２．５％、約２．２％、約７．５％、約３．２％であり、解析結果は理論値とよく一致していることがわかる。

【0082】

図１１は、等方的に２回ずつくりこみを行った場合の解析結果を示す。図１２、図１３、図１４は、それぞれｘ方向、ｙ方向、ｚ方向に２回のくりこみを行い、他の方向にはくりこみを行っていない場合の解析結果を示す。なお、図１１～図１４では、式（３４）のくりこみ変換則を用いている。図１１、図１２、図１３、図１４に示した解析結果の理論値に対する平均誤差は、それぞれ約１．９％、約４．２％、約８．４％、約２．７％であり、解析結果は理論値とよく一致していることがわかる。

【0083】

図１５は、等方的に３回ずつくりこみを行った場合の解析結果を示す。図１６、図１７、図１８は、それぞれｘ方向、ｙ方向、ｚ方向に３回のくりこみを行い、他の方向にはくりこみを行っていない場合の解析結果を示す。なお、図１５～図１８では、式（３４）の変換則を用いている。図１５、図１６、図１７、図１８に示した解析結果の理論値に対する平均誤差は、それぞれ約２．７％、約４．４％、約１０．０％、約７．２％であり、解析結果は理論値とよく一致していることがわかる。

【0084】

図１９～図２２は、式（３３）に示したくりこみ変換則を用いて解析を行った結果を示すグラフである。図１９は、等方的に３回ずつくりこみを行った場合の解析結果を示す。図２０、図２１、図２２は、それぞれｘ方向、ｙ方向、ｚ方向に３回のくりこみを行い、他の方向にはくりこみを行っていない場合の解析結果を示す。図１９、図２０、図２１、図２２に示した解析結果の理論値に対する平均誤差は、それぞれ約４．１％、約３３．６％、約８６．９％、約７１．１％である。なお、図１５と図１９とでは、くりこみ条件は同一であるが、粒子３１の初期の配置、式（３）のランダム力Ｒｉを発生させるときの乱数の違い等の影響により、解析結果に差が生じている。なお、長時間の時間アンサンブルをとれば、両者の違いは小さくなっていくと考えられる。

【0085】

等方的にくりこみを行う場合には、式（３３）の変換則を用いても、解析結果は理論値とよく一致していることがわかる。ところが、異方的にくりこみを行う場合には、式（３３）の変換則を適用すると、解析結果が理論値から大きくずれてしまう。異方的なくりこみを行う場合には、式（３３）の変換則を適用できないことがわかる。

【0086】

図６に示した解析結果から、くりこみを行わない場合に、図１Ａ～図４を参照して説明した実施例によるシミュレーション方法により、壁面４０に接する流体の流れを精度よく解析できることが確認された。図７～図１８に示した解析結果から、等方的または異方的にくりこみを行う場合に、式（３４）の変換則を用いることにより、壁面４０に接する流体の流れを精度よく解析できることが確認された。図１９に示した解析結果から、等方的にくりこみを行う場合に、式（３３）の変換則を用いることにより、壁面４０に接する流体の流れを精度よく解析できることが確認された。

【0087】

次に、図２３Ａ及び図２３Ｂを参照して、さらに他の実施例によるシミュレーション方法及びシミュレーション装置について説明する。

【0088】

図５に示した実施例では、平行な２枚の壁面４０の間を一方向に流れる流体を例にとって、くりこみ手法を説明した。式（３４）に示した変換則は、回転する流れについても適用可能である。

【0089】

図２３Ａ及び図２３Ｂは、本実施例によるシミュレーション方法の解析対象となる解析モデルの断面図である。図２３Ａは、図２３Ｂの一点鎖線２３Ａ－２３Ａにおける断面図であり、図２３Ｂは、図２３Ａの一点鎖線２３Ｂ－２３Ｂにおける断面図である。

【0090】

直径の異なる２つの円筒状の壁面４０Ａ、４０Ｂが同心円状に配置されている。壁面４０Ａ、４０Ｂの間に、流体が周回する流路が形成される。円筒状の壁面４０Ａ、４０Ｂの中心軸に平行な方向をｙ方向とするｘｙｚ直交座標系を定義する。ｘ方向とｚ方向とは等価であるため、くりこみを行う際にｘ方向のくりこみ回数とｚ方向のくりこみ回数とは同一に設定される。このとき、ｘ方向のくりこみ因子λ_ｘとｚ方向のくりこみ因子λ_ｚとは等しい。ｘ方向及びｚ方向のくりこみ因子をλ_ｘｚと標記すると、式（３４）は以下のように変形される。

【数36】

【0091】

式（３６）の変換則を用いることにより、回転流体の解析を行うことが可能である。この変換則は、例えば、外側の壁面４０Ｂをケーシング、内側の壁面４０Ａを攪拌翼として、攪拌槽内の流体の流れの解析に応用することが可能である。

【0092】

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

【符号の説明】

【0093】

３０ 解析空間
３１ 粒子
３１Ｖ 仮想的な粒子
３２、３３ くりこまれた粒子
４０、４０Ａ、４０Ｂ 壁面
５０ 入力部
５１ 処理部
５２ 出力部
５３ 記憶部

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】

【図12】

【図13】

【図14】

【図15】

【図16】

【図17】

【図18】

【図19】

【図20】

【図21】

【図22】

【図23】

【手続補正書】

【提出日】2023-04-05

【手続補正1】

【補正対象書類名】特許請求の範囲

【補正対象項目名】請求項３

【補正方法】変更

【補正の内容】

【請求項3】

壁面に沿って流れる流体の流れを解析するシミュレーション装置であって、
シミュレーション条件が入力される入力部と、
前記入力部に入力された前記シミュレーション条件に基づいて、前記流体の流れを解析する処理部と、
前記処理部による解析結果が出力される出力部と
を備え、
前記処理部は、前記入力部に入力された前記シミュレーション条件に基づいて前記流体を複数の粒子で表し、
前記複数の粒子のそれぞれについて、前記複数の粒子の運動を支配する運動方程式を解くことにより、前記複数の粒子の位置及び速度を時間発展させ、
前記運動方程式を解く際に、
前記複数の粒子のうち前記壁面までの距離が、前記シミュレーション条件で設定された第１距離以下の粒子について、前記シミュレーション条件で設定された粒子間相互作用及び粒子壁面間相互作用による力、前記壁面から受ける減衰力、及び前記壁面の温度に応じたランダム力を作用させて粒子の位置及び速度を時間発展させるシミュレーション装置。

【手続補正2】

【補正対象書類名】特許請求の範囲

【補正対象項目名】請求項５

【補正方法】変更

【補正の内容】

【請求項5】

壁面に沿って流れる流体の流れを解析する手順をコンピュータに実行させるプログラムであって、
シミュレーション条件を取得する手順と、
取得された前記シミュレーション条件に基づいて、前記流体の流れを解析する手順と
を実行させ、
前記流体の流れを解析する手順は、
取得された前記シミュレーション条件に基づいて前記流体を複数の粒子で表す手順と、
前記複数の粒子のそれぞれについて、前記複数の粒子の運動を支配する運動方程式を解くことにより、前記複数の粒子の位置及び速度を時間発展させる手順と
を含み、
前記運動方程式を解く際に、
前記複数の粒子のうち前記壁面までの距離が、前記シミュレーション条件で設定される第１距離以下の粒子について、前記シミュレーション条件で設定された粒子間相互作用及び粒子壁面間相互作用による力、前記壁面から受ける減衰力、及び前記壁面の温度に応じたランダム力を作用させて粒子の位置及び速度を時間発展させるプログラム。

【手続補正3】

【補正対象書類名】明細書

【補正対象項目名】００１０

【補正方法】変更

【補正の内容】

【0010】

【図1】図１Ａ及び図１Ｂは、一実施例によるシミュレーション方法で解析対象となる解析モデルの一例を模式的に示す断面図である。

【図2】図２Ａは、壁面からの距離が第１距離より遠い位置の粒子に作用する力を示す模式図であり、図２Ｂは、壁面からの距離が第１距離以下の位置にある粒子に作用する力を示す模式図である。

【図3】図３は、本実施例によるシミュレーション装置のブロック図である。

【図4】図４は、本実施例によるシミュレーション方法の手順を示すフローチャートである。

【図5】図５Ａは、解析対象の流体を複数の粒子で表した解析モデルの模式図であり、図５Ｂは、図５Ａに示した粒子系に対して、等方的なくりこみを行った粒子系を示す模式図であり、図５Ｃは、図５Ｂに示した粒子系に対して、ｙ方向にはくりこみを行わず、ｘ方向及びｚ方向にさらにくりこみを行った粒子系を示す模式図である。

【図6】図６は、くりこみを行わない解析モデルを用いた解析結果から計算したｘ方向の流速のｙ方向の分布を示すグラフである。

【図7】図７は、ｘ、ｙ、ｚ方向に１回のくりこみを行った解析モデルを用いた解析結果から計算したｘ方向の流速のｙ方向の分布を示すグラフである。

【図8】図８は、ｘ方向に１回のくりこみを行った解析モデルを用いた解析結果から計算したｘ方向の流速のｙ方向の分布を示すグラフである。

【図9】図９は、ｙ方向に１回のくりこみを行った解析モデルを用いた解析結果から計算したｘ方向の流速のｙ方向の分布を示すグラフである。

【図10】図１０は、ｚ方向に１回のくりこみを行った解析モデルを用いた解析結果から計算したｘ方向の流速のｙ方向の分布を示すグラフである。

【図11】図１１は、ｘ、ｙ、ｚ方向に２回のくりこみを行った解析モデルを用いた解析結果から計算したｘ方向の流速のｙ方向の分布を示すグラフである。

【図12】図１２は、ｘ方向に２回のくりこみを行った解析モデルを用いた解析結果から計算したｘ方向の流速のｙ方向の分布を示すグラフである。

【図13】図１３は、ｙ方向に２回のくりこみを行った解析モデルを用いた解析結果から計算したｘ方向の流速のｙ方向の分布を示すグラフである。

【図14】図１４は、ｚ方向に２回のくりこみを行った解析モデルを用いた解析結果から計算したｘ方向の流速のｙ方向の分布を示すグラフである。

【図15】図１５は、ｘ、ｙ、ｚ方向に３回のくりこみを行った解析モデルを用いた解析結果から計算したｘ方向の流速のｙ方向の分布を示すグラフである。

【図16】図１６は、ｘ方向に３回のくりこみを行った解析モデルを用いた解析結果から計算したｘ方向の流速のｙ方向の分布を示すグラフである。

【図17】図１７は、ｙ方向に３回のくりこみを行った解析モデルを用いた解析結果から計算したｘ方向の流速のｙ方向の分布を示すグラフである。

【図18】図１８は、ｚ方向に３回のくりこみを行った解析モデルを用いた解析結果から計算したｘ方向の流速のｙ方向の分布を示すグラフである。

【図19】図１９は、ｘ、ｙ、ｚ方向に３回のくりこみを行った解析モデルを用いた解析結果から計算したｘ方向の流速のｙ方向の分布を示すグラフである。

【図20】図２０は、ｘ方向に３回のくりこみを行った解析モデルを用い、比較例による方法で解析を行った結果から計算したｘ方向の流速のｙ方向の分布を示すグラフである。

【図21】図２１は、ｙ方向に３回のくりこみを行った解析モデルを用い、比較例による方法で解析を行った結果から計算したｘ方向の流速のｙ方向の分布を示すグラフである。

【図22】図２２は、ｚ方向に３回のくりこみを行った解析モデルを用い、比較例による方法で解析を行った結果から計算したｘ方向の流速のｙ方向の分布を示すグラフである。

【図23】図２３Ａ及び図２３Ｂは、他の実施例によるシミュレーション方法の解析対象となる解析モデルの断面図である。

【手続補正4】

【補正対象書類名】明細書

【補正対象項目名】００２２

【補正方法】変更

【補正の内容】

【0022】

減衰係数γは、例えば以下の式で表される。

【数5】

ここで、ζは減衰比であり、ｋはバネ定数である。減衰比ζとして、例えば０．７０７を用いることができる。バネ定数ｋは、例えば、レナードジョーンズポテンシャル（式（１））を鞍点近傍でテイラー展開したときの２次の項の係数から求めることができる。この場合、バネ定数ｋは、以下の式で計算することができる。

【数6】

減衰係数γとして、例えば５．２３×１０ ^－１３ ｋｇ／ｓを用いることができる。

【手続補正5】

【補正対象書類名】明細書

【補正対象項目名】００４４

【補正方法】変更

【補正の内容】

【0044】

従って、運動エネルギＥ_ｋは、以下の式で表される。

【数12】

くりこみ因子λのくりこみを行うと、式（１２）の右辺のＮ及びβが共に１／λ^３になるため、運動エネルギＥ_Ｋは、くりこみの前後で不変である。

【手続補正6】

【補正対象書類名】明細書

【補正対象項目名】００４８

【補正方法】変更

【補正の内容】

【0048】

従って、相互作用エネルギＥ_ｉｎｔは、以下の式で近似することができる。

【数16】

【手続補正7】

【補正対象書類名】明細書

【補正対象項目名】００６９

【補正方法】変更

【補正の内容】

【0069】

ｘ方向、ｙ方向、及びｚ方向へのくりこみ因子を、それぞれλ_ｘ、λ_ｙ、及びλ_ｚと標記する。このときのくりこみ変換則は、以下の式で表される。

【数28】

【手続補正8】

【補正対象書類名】明細書

【補正対象項目名】００７７

【補正方法】変更

【補正の内容】

【0077】

次に、図６～図２２を参照して、実際にシミュレーションを行った結果について説明する。図１Ａ、図１Ｂに示したポアズイユ流れの解析モデルについて、くりこみを行わない方法、等方的なくりこみを行う方法、異方的なくりこみを行う方法で解析を行った。ｘ方向とｚ方向とに、周期境界条件を適用した。ｘ方向の寸法Ｌ _ｘ 及びｚ方向の寸法Ｌ _ｚ をともに７．２０ｎｍとし、ｙ方向の寸法Ｌ _ｙ を６．２９ｎｍとした。