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

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

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

(11)【公開番号】P2024175843

(43)【公開日】2024-12-19

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

(51)【国際特許分類】

   G16Z 99/00 20190101AFI20241212BHJP

   G01N 11/04 20060101ALN20241212BHJP

【ＦＩ】

G16Z99/00

G01N11/04 Z

【審査請求】未請求

【請求項の数】6

【出願形態】ＯＬ

(21)【出願番号】P 2023093889

(22)【出願日】2023-06-07

(71)【出願人】

【識別番号】000002107

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

(74)【代理人】

【識別番号】100105887

【弁理士】

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

(72)【発明者】

【氏名】廣田 真人

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

5L049

【Ｆターム（参考）】

5L049AA24

(57)【要約】

【課題】解析の前処理として、解析空間内の領域ごとの空間解像度の設定を行う必要が無く、かつ計算負荷を低減させることが可能なシミュレーション装置を提供する。
【解決手段】流体を複数の計算粒子で表し、複数の計算粒子のそれぞれに物理量を付与し、流路をモデル化した解析空間に複数の計算粒子を配置し、複数の計算粒子のそれぞれに付与された物理量及び複数の計算粒子の位置を、流体の支配方程式を解いて時間発展させる。複数の計算粒子のそれぞれに付与された物理量及び複数の計算粒子の位置を時間発展させる処理において、複数の計算粒子のそれぞれの位置における流路の特徴的な寸法と、複数の計算粒子のそれぞれの大きさとを比較し、比較結果に応じて複数の計算粒子のそれぞれを分割して、より小さな計算粒子を生成する処理、または複数の計算粒子のいくつかを合体させてより大きな計算粒子を生成する処理を行う。
【選択図】図２

【特許請求の範囲】

【請求項1】

流路の内部の流体の流れを、粒子法を用いて解析するシミュレーション装置であって、
シミュレーション条件が入力される入力部と、
前記入力部に入力されたシミュレーション条件に基づいて演算を行う処理部と
を備え、
前記処理部は、
流体を複数の計算粒子で表し、前記複数の計算粒子のそれぞれに物理量を付与し、前記流路をモデル化した解析空間に前記複数の計算粒子を配置し、前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を、流体の支配方程式を解いて時間発展させる機能と、
前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を時間発展させる処理において、前記複数の計算粒子のそれぞれの位置における前記流路の特徴的な寸法と、前記複数の計算粒子のそれぞれの大きさとを比較し、比較結果に応じて前記複数の計算粒子のそれぞれを分割して、より小さな計算粒子を生成する処理、または前記複数の計算粒子のいくつかを合体させてより大きな計算粒子を生成する処理を行う機能と
を有するシミュレーション装置。

【請求項2】

前記複数の計算粒子のそれぞれの位置における前記流路の特徴的な寸法は、前記複数の計算粒子のそれぞれの位置における前記流路の空間的な広がりの程度に基づいて求められる請求項１に記載のシミュレーション装置。

【請求項3】

前記解析空間に張られた格子の複数の格子点のそれぞれに対して、前記流路の壁面までの最短距離の値が割り当てられており、
前記処理部は、前記複数の計算粒子のそれぞれの位置の近傍に位置する格子点に割り当てられた前記最短距離の値に基づいて、前記複数の計算粒子のそれぞれの位置における前記特徴的な寸法を求める機能を、さらに有する請求項１または２に記載のシミュレーション装置。

【請求項4】

前記処理部は、
前記複数の計算粒子のそれぞれの位置から、前記複数の格子点に割り当てられた前記最短距離の値の勾配を表す方向に、前記最短距離の値の極大値を探索し、探索によって得られた前記極大値に基づいて、前記特徴的な寸法を求める請求項３に記載のシミュレーション装置。

【請求項5】

流路の内部の流体の流れを、粒子法を用いて解析するシミュレーション方法であって、
シミュレーション条件を与えられたコンピュータが、
与えられたシミュレーション条件に基づいて、流体を複数の計算粒子で表し、前記複数の計算粒子のそれぞれに物理量を付与し、前記流路をモデル化した解析空間に前記複数の計算粒子を配置し、前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を、流体の支配方程式を解いて時間発展させ、
前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を時間発展させる処理において、前記複数の計算粒子のそれぞれの位置における前記流路の特徴的な寸法と、前記複数の計算粒子のそれぞれの大きさとを比較し、比較結果に応じて前記複数の計算粒子のそれぞれを分割して、より小さな計算粒子を生成する処理、または前記複数の計算粒子のいくつかを合体させてより大きな計算粒子を生成する処理を行うシミュレーション方法。

【請求項6】

流路の内部の流体の流れを、粒子法を用いて解析する機能をコンピュータに実現させるプログラムであって、
与えられたシミュレーション条件に基づいて、流体を複数の計算粒子で表し、前記複数の計算粒子のそれぞれに物理量を付与し、前記流路をモデル化した解析空間に前記複数の計算粒子を配置し、前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を、流体の支配方程式を解いて時間発展させる機能と、
前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を時間発展させる処理において、前記複数の計算粒子のそれぞれの位置における前記流路の特徴的な寸法と、前記複数の計算粒子のそれぞれの大きさとを比較し、比較結果に応じて前記複数の計算粒子のそれぞれを分割して、より小さな計算粒子を生成する処理、または前記複数の計算粒子のいくつかを合体させてより大きな計算粒子を生成する処理を行う機能と
をコンピュータに実現させるプログラム。

【発明の詳細な説明】

【技術分野】

【0001】

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

【背景技術】

【0002】

流体解析技術として、連続な場を仮想的な粒子（以下、計算粒子という。）の集合を用いて離散化して支配方程式を解くことにより、流れ場の時間発展を解析する粒子法が知られている。この解析において、空間解像度が計算粒子の直径に依存する。従来技術においては、解析の前処理の段階で計算粒子の直径を設定し、解析空間に同じ直径を有する複数の計算粒子を配置して解析を行う。

【0003】

与えられた形状モデル及び計算粒子の直径の最小値に基づいて、解析に用いる計算粒子の直径を設定したり形状モデルを修正したりする方法が知られている（特許文献１）。十分な空間解像度を得るために、最も流路の狭い領域に合わせて計算粒子の直径を設定する必要がある。解析空間に単一の直径を有する計算粒子を配置する場合、流路の狭い領域が解析空間全体に対してわずかであったとしても、解析空間全体において流路の狭い領域に合わせて計算粒子の直径を設定しなければならない。このとき、流路が広い領域では空間解像度が過剰になり、解析負荷が増加してしまう。

【0004】

解析空間を複数の領域に分割し、分割された領域ごとに異なる粒子径の計算粒子を配置する手法が公知である（特許文献２）。また、解析空間を壁からの距離に応じて複数の領域に分割し、分割された領域ごとに異なる粒子径の計算粒子を配置する手法が公知である（特許文献３）。

【先行技術文献】

【特許文献】

【0005】

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

【特許文献2】特開２００７－１４０９９３号公報

【特許文献3】特開２０２２－１０８６５８号公報

【発明の概要】

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

【0006】

従来の手法では、高い空間解像度を要する領域とそうでない領域とをオペレータが事前に把握する必要があり、解析の前処理における工数が増える。また、オペレータの技量によっては、空間解像度の不適切な設定がなされる可能性がある。本発明の目的は、解析の前処理として、解析空間内の領域ごとの空間解像度の設定を行う必要が無く、かつ計算負荷を低減させることが可能なシミュレーション装置、シミュレーション方法、及びプログラムを提供することである。

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

【0007】

本発明の一観点によると、
流路の内部の流体の流れを、粒子法を用いて解析するシミュレーション装置であって、
シミュレーション条件が入力される入力部と、
前記入力部に入力されたシミュレーション条件に基づいて演算を行う処理部と
を備え、
前記処理部は、
流体を複数の計算粒子で表し、前記複数の計算粒子のそれぞれに物理量を付与し、前記流路をモデル化した解析空間に前記複数の計算粒子を配置し、前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を、流体の支配方程式を解いて時間発展させる機能と、
前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を時間発展させる処理において、前記複数の計算粒子のそれぞれの位置における前記流路の特徴的な寸法と、前記複数の計算粒子のそれぞれの大きさとを比較し、比較結果に応じて前記複数の計算粒子のそれぞれを分割して、より小さな計算粒子を生成する処理、または前記複数の計算粒子のいくつかを合体させてより大きな計算粒子を生成する処理を行う機能と
を有するシミュレーション装置が提供される。

【0008】

本発明の他の観点によると、
流路の内部の流体の流れを、粒子法を用いて解析するシミュレーション方法であって、
シミュレーション条件を与えられたコンピュータが、
与えられたシミュレーション条件に基づいて、流体を複数の計算粒子で表し、前記複数の計算粒子のそれぞれに物理量を付与し、前記流路をモデル化した解析空間に前記複数の計算粒子を配置し、前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を、流体の支配方程式を解いて時間発展させ、
前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を時間発展させる処理において、前記複数の計算粒子のそれぞれの位置における前記流路の特徴的な寸法と、前記複数の計算粒子のそれぞれの大きさとを比較し、比較結果に応じて前記複数の計算粒子のそれぞれを分割して、より小さな計算粒子を生成する処理、または前記複数の計算粒子のいくつかを合体させてより大きな計算粒子を生成する処理を行うシミュレーション方法が提供される。

【0009】

本発明の他の観点によると、
流路の内部の流体の流れを、粒子法を用いて解析する機能をコンピュータに実現させるプログラムであって、
与えられたシミュレーション条件に基づいて、流体を複数の計算粒子で表し、前記複数の計算粒子のそれぞれに物理量を付与し、前記流路をモデル化した解析空間に前記複数の計算粒子を配置し、前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を、流体の支配方程式を解いて時間発展させる機能と、
前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を時間発展させる処理において、前記複数の計算粒子のそれぞれの位置における前記流路の特徴的な寸法と、前記複数の計算粒子のそれぞれの大きさとを比較し、比較結果に応じて前記複数の計算粒子のそれぞれを分割して、より小さな計算粒子を生成する処理、または前記複数の計算粒子のいくつかを合体させてより大きな計算粒子を生成する処理を行う機能と
をコンピュータに実現させるプログラムが提供される。

【発明の効果】

【0010】

流路の特徴的な寸法に応じて計算粒子の合体または分割が行われるため、解析の前処理として、解析空間内の領域ごとの空間解像度の設定を行う必要が無い。さらに、計算粒子が合体されると、計算負荷を低減させることができる。

【図面の簡単な説明】

【0011】

【図1】図１Ａは、解析対象の流体が流れる流路の形状の一例を示す断面図であり、図１Ｂは、符号付距離関数の絶対値を濃淡で表した図であり、図１Ｃは、流路の任意の位置における特徴的な寸法を求める方法を説明するための模式図であり、図１Ｄは、算出された特徴的な寸法の分布を、濃淡で表した図である。

【図2】図２は、流れの解析中における計算粒子の合体及び分割の様子の一例を示す模式図である。

【図3】図３Ａは、分割前の計算粒子の模式図であり、図３Ｂは、１つの計算粒子を分割して生成された２個の計算粒子の模式図である。

【図4】図４Ａは、１つの計算粒子を分割して生成された４個の計算粒子の模式図であり、図４Ｂは、１つの計算粒子を分割して生成された８個の計算粒子の模式図である。

【図5】図５Ａは、複数の計算粒子の分布の一例を示す模式図であり、図５Ｂは、２つの計算粒子を合体させた後の計算粒子の分布の一例を示す模式図である。

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

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

【図8】図８Ａは、計算開始時点における計算粒子の分布を示す図であり、図８Ｂは、定常状態における計算粒子の分布を示す図である。

【図9】図９は、定常状態における流路の中心軸上での圧力及び液滴の平均半径の時間平均値の分布の計算結果を、論文に示された実験結果と比較して示すグラフである。

【図10】図１０Ａ及び図１０Ｂは、変形例によるシミュレーション方法で解析対象となる流路の、流れ方向に直交する断面を示す図である。

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

【0012】

［粒子法の概要］
本発明の一実施例について説明する前に、一般的なＳＰＨ（Smoothed Particle Hydrodynamics）法について簡単に説明する。ＳＰＨ法は、連続体を複数の計算粒子の集合を用いて離散化し、流れ場の時間発展を計算するラグランジュ型の数値計算手法である。ＳＰＨ法では、粒子１個の質量分布として、連続かつ微分可能なカーネル関数を用いるとともに、カーネル関数の微分係数の重ね合わせを空間微分とする離散化モデルを用いる。

【0013】

ＳＰＨ法では、計算粒子の質量ｍはカーネル関数Ｗによってその周囲に分配され、連続的に分布しており、密度場は、その重ね合わせとして与えられる。すなわち、ｉ番目の計算粒子の位置ｘ_ｉにおける密度ρ_ｉは、以下の式で表される。

【数1】

本明細書において、下付きの添え字ｉ及びｊは、計算粒子のそれぞれに付与された通し番号を表し、下付きの添え字ｉ、ｊが付された物理量は、それぞれｉ番目及びｊ番目の計算粒子に付された物理量を意味する。

【0014】

カーネル関数Ｗとして、例えば、以下の関数を用いることができる。以下の関数は、Ｃ２Ｗｅｎｄｌａｎｄカーネルと呼ばれる。

【数2】

【0015】

ｈは、平滑化長さであり、カーネル幅ともいわれる。カーネル幅ｈは、複数の計算粒子のそれぞれを同一の質量及び物性密度の球体に置き換えた等価球の直径に等しい。本明細書において、カーネル幅ｈを計算粒子の直径という場合がある。

【0016】

任意の位置ｘ_ｉにおける任意のスカラー量Ａ及びベクトル量Ａの値は、カーネル関数の重ね合わせによって以下の式で表すことができる。

【数3】

【0017】

任意のスカラー量Ａの勾配は、カーネル関数の微分係数の重ね合わせによって以下の式で表される。

【数4】

【数5】

【0018】

ここで、∇Ｗ（ｘ_ｉｊ，ｈ）はカーネル関数の導関数であり、カーネル関数としてＣ２Ｗｅｎｄｌａｎｄカーネルを用いる場合、以下の式で表される。

【数6】

なお、式（２）の関数Ｗ（ｑ）は１変数の関数であるため、その勾配∇Ｗはスカラー量である。

【0019】

任意のベクトル量Ａに対する発散は、式（４）と同様に、以下の式で表すことができる。

【数7】

【0020】

任意のスカラー量Ａに対する２次導関数は、例えば以下の公式を用いて算出することができる。

【数8】

【0021】

［湿り蒸気流れの支配方程式］
非粘性、断熱、均質流を仮定した相変化を伴う圧縮性流れの支配方程式は、以下のように表される。

【数9】

【数10】

【数11】

【数12】

【数13】

ここで、ｔは時間、ρは二相流体の密度、ｖは速度ベクトル、ｐは圧力、ｕは内部エネルギを表す。Ψは液相の質量分率、Γは液相の質量生成率、ｎは液滴の数密度、Ｉは凝縮核生成率を表す。

【0022】

上述の式（９）～（１３）を、ＳＰＨ法を用いて離散化すると、以下のように表すことができる。

【数14】

【数15】

【数16】

【数17】

【数18】

【数19】

ここで、下付き文字のｉ及びｊは、計算粒子に付された番号を表す。Π_ｉｊは人工粘性項であり、以下の式で表される。

【数20】

【数21】

ここで、εは、ｘ_ｉｊ＝０のときの特異性を防ぐための定数であり、例えばε＝０．０１とすることができる。ａは、人工粘性の強さを調節するパラメータであり、例えばａ＝１とすることができる。

【0023】

［流路の特徴的な寸法の算出］
本実施例では、流体が流れる流路の特徴的な寸法に基づいて、計算途中に計算粒子を分割したり、合体させたりする。次に、図１Ａ～図１Ｄを参照して流路の特徴的な寸法について説明する。

【0024】

図１Ａは、解析対象の流体が流れる流路１０の形状の一例を示す断面図である。物体Ωによって流路１０の形状及び大きさが決まる。物体Ωの形状及び大きさは、例えばＣＡＤデータで与えられる。物体Ωによって画定される流路１０の表面を∂Ωと標記する。図１Ａにおいては、二次元の解析空間内に物体Ωが配置され、二次元の流路１０が定義されている。

【0025】

図１Ａにおいて、流路１０の左端から流体が流入し、右端から流出する。流路１０の流入端の位置の流路幅をＷｉと標記し、流出端の位置の流路幅をＷｏと標記し、流路１０の長さをＬと標記する。流路は、流入端の近傍に流路幅一定の部分を有し、下流に向って流路幅が徐々に狭くなってスロート部を形成し、その後徐々に広くなる。流路幅が徐々に狭くなる部分の長さは、流路幅が徐々に広くなる部分の長さより短い。流出端の近傍では、流路幅が一定になる。流出端の位置における流路幅Ｗｏは、流入端の位置における流路幅Ｗｉより小さい。スロート部の流路幅をＷｔと標記する。

【0026】

粒子法を用いた流体解析において、物体周りの流れを解析する際に、物体表面における境界条件を与えるために、複数の計算粒子のそれぞれの位置から物体表面までの距離を算出する必要がある。

【0027】

解析空間の任意の位置における物体表面からの距離を算出するために、例えば符号付距離関数（Signed Distance Function）を用いることができる。解析空間の任意の位置をｘ、位置ｘにおける符号付距離関数の値をｆ（ｘ）と標記すると、ｆ（ｘ）は以下のように定義される。

【数22】

ここで、ｄ（ｘ，∂Ω）は、位置ｘにおける表面∂Ωからの最短距離を表す。すなわち、符号付距離関数の絶対値は、物体表面からの最短距離であり、符号付距離関数の符号は、物体の外部で正、物体の内部で負である。符号付距離関数は、計算力学や可視化情報学の分野で広く用いられている。

【0028】

図１Ｂは、符号付距離関数の絶対値を濃淡で表した図である。物体Ωの表面∂Ωにおいて符号付距離関数の値がゼロになり、表面∂Ωから遠ざかるにしたがって符号付距離関数の絶対値が大きくなる（図１Ｂにおいて濃くなる）。流路１０の内部（物体Ωの外部）では、符号付距離関数の値が正であり、流路１０の外部（物体Ωの内部）では、符号付距離関数の値が負である。図１Ｂに示した破線は、符号付距離関数の値の等値線を表す。

【0029】

図１Ｂに示した細い実線は、解析空間に張られた直交等間隔格子を表す。直交等間隔格子の複数の格子点のそれぞれの位置において、符号付距離関数の値が求められている。なお、実際の解析で用いる直交等間隔格子の格子間隔は、図１Ｂに示した直交等間隔格子の格子間隔より狭い。

【0030】

図１Ｃは、流路１０の任意の位置における特徴的な寸法を求める方法を説明するための模式図である。特徴的な寸法は、解析空間に張られた直交等間隔格子の格子点のそれぞれの位置において計算される。まず、格子点のそれぞれの位置において、符号付距離関数の値の勾配ベクトルを計算する。図１Ｃにおいて、いくつかの格子点のそれぞれの位置における勾配ベクトルを矢印で表している。解析空間内の任意の位置における勾配ベクトルは、符号付距離関数の値の等値線に対して垂直であり、符号付距離関数の値が大きい方を向く。

【0031】

次に、求められた勾配ベクトルの方向に沿って、符号付距離関数の値の極大値を探索する。例えば、極大値を探索するときの勾配ベクトルに沿う直線上の複数の計算位置の刻み幅は、直交等間隔格子の格子間隔と等しくする。図１Ｃに示した勾配ベクトルを示す矢印の始点が、特徴的な寸法を求める位置に相当し、終点が、符号付距離関数の値が極大値を取る位置に相当する。探索された極大値の２倍の値を、格子点のそれぞれの位置における特徴的な寸法として採用する。

【0032】

なお、符号付距離関数の値が負の格子点（物体Ωの内部の格子点）については、上述の計算は行わず、特徴的な寸法としてゼロを割り当てる。

【0033】

図１Ｄは、算出された特徴的な寸法の分布を、濃淡で表した図である。図１Ｄにおいて、領域Ａ１内の格子点における特徴的な寸法が相対的に大きく、領域Ａ３内の格子点における特徴的な寸法が相対的に小さい。領域Ａ２、Ａ４内の格子点における特徴的な寸法は、領域Ａ１内の大きさと領域Ａ４内の大きさとの中間である。すなわち、二次元の流路１０において、流路幅が相対的に広い領域の格子点における特徴的な寸法が、相対的に大きく、流路幅が相対的に狭い領域の格子点における特徴的な寸法が、相対的に小さい。すなわち、特徴的な寸法は、流路幅の値を反映したものということができる。

【0034】

［計算粒子の合体及び分割］
次に、計算粒子を合体させたり分割したりする条件について説明する。
複数の計算粒子のそれぞれの位置における流路の特徴的な寸法を計算する。例えば、解析空間に張られた直交等間隔格子の格子点のそれぞれに割り当てられた流路の特徴的な寸法の値を用いて線形補間演算を行うことにより、複数の計算粒子のそれぞれの位置における流路の特徴的な寸法を計算することができる。このとき、格子点に割り当てられた符号付距離関数の値が負である場合は、その格子点を補間演算の対象から除外する。

【0035】

着目する計算粒子に付与された通し番号をｉと標記し、ｉ番目の計算粒子の位置における流路の特徴的な寸法の値をｗ_ｉと標記し、ｉ番目の計算粒子のカーネル幅（直径）をｈ_ｉと標記する。

【0036】

以下の式（２３）が満たされる場合、空間解像度が不足していると判定し、計算粒子を分割する。

【数23】

ここで、α_ｓは調整パラメータであり、例えばα_ｓとして０．２を用いることができる。すなわち、計算粒子の直径ｈ_ｉが、その位置における流路の特徴的な寸法ｗ_ｉにある定数α_ｓを乗じた値より大きい場合、空間解像度が不足していると判定される。

【0037】

以下の式（２４）が満たされる場合、空間解像度が過剰であると判定し、計算粒子を合体させる。

【数24】

ここで、α_ｍは調整パラメータであり、例えばα_ｍとして０．１を用いることができる。すなわち、計算粒子の直径ｈ_ｉが、その位置における流路の特徴的な寸法ｗ_ｉにある定数α_ｍを乗じた値より小さい場合、空間解像度が過剰であると判定される。なお、定数α_ｍを定数α_ｓより小さくすることが好ましい。

【0038】

図２は、流れ場の解析中における計算粒子の合体及び分割の様子の一例を示す模式図である。図１Ｄを参照して説明したように、領域Ａ１内において特徴的な寸法が相対的に大きく、領域Ａ３内において特徴的な寸法が相対的に小さく、領域Ａ２、Ａ４内において特徴的な寸法が、その中間の大きさである。

【0039】

流路１０の流入端（図２において左端）から流出端（図２において右端）に向って計算粒子２０が移動する。計算粒子２０が領域Ａ１から領域Ａ２に進入すると、特徴的な寸法の値ｗ_ｉが小さくなる。これにより式（２３）が満たされると、計算粒子２０が複数の計算粒子２０に分割されて、計算粒子２０のそれぞれの直径が小さくなる。同様に、計算粒子２０が領域Ａ２から領域Ａ３に進入するときも、特徴的な寸法の値ｗ_ｉが小さくなる。これにより式（２３）が満たされると、計算粒子２０がさらに分割されて、計算粒子２０の直径がさらに小さくなる。

【0040】

計算粒子２０が領域Ａ３から領域Ａ４に進入すると、特徴的な寸法の値ｗ_ｉが大きくなる。これにより式（２４）が満たされると、複数の計算粒子２０が１つの計算粒子２０に合体されて、計算粒子２０の直径が大きくなる。

【0041】

次に、図３Ａ～図４Ｂを参照して、計算粒子を分割する方法について説明する。
図３Ａは、分割前の計算粒子２０ａの模式図である。図３Ｂは、１つの計算粒子２０ａ（図３Ａ）を分割して生成された２個の計算粒子２０ａ１、２０ａ２の模式図である。図４Ａは、１つの計算粒子２０ａ（図３Ａ）を分割して生成された４個の計算粒子２０ａ１～２０ａ４の模式図である。図４Ｂは、１つの計算粒子２０ａ（図３Ａ）を分割して生成された８個の計算粒子２０ａ１～２０ａ８の模式図である。

【0042】

まず、図３Ｂに示すように、２個の計算粒子２０ａ１、２０ａ２に分割する場合について説明する。分割前の計算粒子２０ａの位置をｘ_１、カーネル幅をｈ_１、質量をｍ_１、体積をＶ_１と標記する。分割後の２つの計算粒子２０ａ１、２０ａ２の各々のカーネル幅をｈ_２、質量をｍ_２、体積をＶ_２と標記する。分割後の２つの計算粒子２０ａ１、２０ａ２の位置を、それぞれｘ_２１、ｘ_２２と標記する。ここで、Ｖ_１＝ｈ_１ ^３、Ｖ_２＝ｈ_２ ^３である。

【0043】

分割後の計算粒子２０ａ１、２０ａ２のそれぞれの質量ｍ_２、体積Ｖ_２、カーネル幅ｈ_２は、以下の式で計算することができる。

【数25】

【0044】

分割後の２つの計算粒子２０ａ１、２０ａ２の位置ｘ_２１、ｘ_２２は、例えば以下の式で与えることができる。

【数26】

ここで、ベクトルｉ、ｊ、ｋは、それぞれ解析空間に定義されたｘｙｚ直交座標系のｘ方向、ｙ方向、ｚ方向の単位ベクトルを表す。

【0045】

次に、図４Ａに示すように、４個の計算粒子２０ａ１～２０ａ４に分割する場合について説明する。分割後の４個の計算粒子２０ａ１～２０ａ４の各々のカーネル幅をｈ_２、質量をｍ_２、体積をＶ_２と標記する。分割後の４個の計算粒子２０ａ１～２０ａ４の位置を、それぞれｘ_２１～ｘ_２４と標記する。

【0046】

分割後の計算粒子２０ａ１～２０ａ４のそれぞれの質量ｍ_２、体積Ｖ_２、カーネル幅ｈ_２は、以下の式で計算することができる。

【数27】

【0047】

分割後の４個の計算粒子２０ａ１～２０ａ４の位置ｘ_２１～ｘ_２４は、例えば以下の式で与えることができる。

【数28】

【0048】

次に、図４Ｂに示すように、８個の計算粒子２０ａ１～２０ａ８に分割する場合について説明する。分割後の８個の計算粒子２０ａ１～２０ａ８の各々のカーネル幅をｈ_２、質量をｍ_２、体積をＶ_２と標記する。分割後の８個の計算粒子２０ａ１～２０ａ８の位置を、それぞれｘ_２１～ｘ_２８と標記する。

【0049】

分割後の計算粒子２０ａ１～２０ａ８のそれぞれの質量ｍ_２、体積Ｖ_２、カーネル幅ｈ_２は、以下の式で計算することができる。

【数29】

【0050】

分割後の８個の計算粒子２０ａ１～２０ａ８の位置ｘ_２１～ｘ_２８は、例えば以下の式で与えることができる。

【数30】

【0051】

分割後のｓ番目の計算粒子２０ｓの任意の物理量Ａ_２ｓは、例えば分割前の計算粒子２０の物理量Ａ_１の勾配を用いて以下の式で与えることができる。

【数31】

ここで、物理量Ａ_２ｓ、Ａ_１は、式（１４）～式（１９）の支配方程式の左辺の物理量に相当する。

【0052】

なお、解析空間が二次元である場合も、同様の手法で分割後の計算粒子２０の位置、及び計算粒子２０の物理量を求めることができる。

【0053】

次に、図５Ａ、図５Ｂを参照して、計算粒子を合体させる方法について説明する。計算粒子を合体させる場合、最近接ペアを構成する相手粒子があるか否かを判定する。「最近接ペア」とは、互いに最近接の関係を有する計算粒子のペアを意味する。

【0054】

図５Ａを参照して、最近接ペアについて具体的に説明する。図５Ａは、複数の計算粒子２０の分布の一例を示す模式図である。図５Ａに示した矢印は、計算粒子２０のそれぞれから最近接の計算粒子２０に向かう。例えば、計算粒子２０ａの最近接粒子は、計算粒子２０ｂであり、計算粒子２０ｂの最近接粒子は計算粒子２０ａである。したがって、２つの計算粒子２０ａ、２０ｂは、最近接ペアを構成する。

【0055】

計算粒子２０ｃの最近接粒子は計算粒子２０ｄである。ところが、計算粒子２０ｄの最近接粒子は、計算粒子２０ｅであり、計算粒子２０ｃではない。したがって、計算粒子２０ｃは、他の計算粒子２０と最近接ペアを構成しない。計算粒子２０ｅの最近接粒子は計算粒子２０ｄである。このため、計算粒子２０ｄと２０ｅとが、最近接ペアを構成している。

【0056】

最近接ペアを構成する相手の計算粒子２０が有る場合、最近接ペアを構成する相手粒子の位置における空間解像度が過剰か否かを判定する。相手粒子の空間解像度が過剰であると判定された場合は、最近接ペアを構成する２つの計算粒子２０を合体させ、２つの計算粒子２０のいずれのカーネル幅より大きいカーネル幅を持つ１つの大きな計算粒子を生成する。すなわち、最近接ペアを構成する２つの計算粒子２０のそれぞれの空間解像度が共に過剰である場合に、２つの計算粒子２０を合体させる。

【0057】

図５Ａにおいて、空間解像度が過剰と判定された計算粒子２０にハッチングを付している。すなわち計算粒子２０ａ、２０ｂ、２０ｃ、２０ｄについて、空間解像度が過剰と判定されたものであり、他の計算粒子２０については、空間解像度が過剰と判定されていない。

【0058】

最近接ペアを構成する計算粒子２０ａと２０ｂとは、共に空間解像度が過剰と判定されたものであるため、計算粒子２０ａと２０ｂとを合体させる。最近接ペアを構成する計算粒子２０ｄと２０ｅとは、一方の計算粒子２０ｅの空間解像度が過剰ではないため、計算粒子２０ｄと２０ｅとは合体させない。

【0059】

図５Ｂは、計算粒子２０ａと２０ｂとを合体させた後の計算粒子２０の分布の一例を示す模式図である。計算粒子２０ａと２０ｂとが合体されて、大きな計算粒子２０ａｂが生成されている。図５Ｂにおいて、合体後の計算粒子２０ａｂに、相対的に淡いハッチングを付している。

【0060】

次に、２つの計算粒子２０を合体させる方法について説明する。
合体前の一方の計算粒子２０ａの位置をｘ_１１、カーネル幅をｈ_１１、質量をｍ_１１、体積をＶ_１１と標記し、他方の計算粒子２０ｂの位置をｘ_１２、カーネル幅をｈ_１２、質量をｍ_１２、体積をＶ_１２と標記する。合体後の計算粒子２０ａｂ（図５Ｂ）の質量ｍ_２、カーネル幅ｈ_２は、以下の式を用いて計算することができる。

【数32】

【0061】

合体後の計算粒子２０ａｂの位置ｘ_２は、合体前の２つの計算粒子２０ａ、２０ｂの位置ｘ_１１、ｘ_１２を、質量に基づく重み付き平均値を用いて、以下の式で決めることができる。

【数33】

【0062】

合体後の計算粒子２０ａｂの任意の物理量Ａ_２も同様に、合体前の計算粒子２０ａ、２０ｂの物理量Ａ_１１、Ａ_１２を用いて以下の式で計算することができる。

【数34】

または、合体前の計算粒子２０ａ、２０ｂの物理量Ａ_１１、Ａ_１２と、これらの物理量の勾配とを用いて、例えば以下の式で計算してもよい。

【数35】

【0063】

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

【0064】

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

【0065】

図７は、本実施例によるシミュレーション方法の手順を示すフローチャートである。図７に示したフローチャートの各ステップは、処理部３１によって実行される。

【0066】

まず、処理部３１が、オペレータによって入力部３０に入力されたシミュレーション条件を取得する（ステップＳ１）。シミュレーション条件には、シミュレーション対象の流体を定義する情報、初期条件、境界条件等が含まれる。流体を定義する情報には、流体の密度、粘度等の物性値が含まれる。初期条件には、複数の計算粒子を解析空間に配置するための情報、計算粒子の初期物性値を指定する方法等が含まれる。境界条件には、解析空間の形状及び大きさを指定する情報、流路を画定する物体の形状を指定する情報等が含まれる。

【0067】

次に、処理部３１は、与えられた物体形状について、直交等間隔格子の格子点のそれぞれの位置における符号付距離関数の値を計算する（ステップＳ２）。さらに、符号付距離関数の値を用いて、格子点のそれぞれの位置における流路１０の特徴的な寸法を計算する（ステップＳ３）。この計算には、図１Ｃを参照して説明した方法を用いることができる。この計算により、図１Ｄに示した特徴的な寸法の分布が得られる。

【0068】

入力されたシミュレーション条件及び符号付距離関数の値に基づいて、解析空間に複数の計算粒子を配置する（ステップＳ４）。具体的には、符号付距離関数の値が正の領域（流路１０の内部）に、複数の計算粒子を配置する。

【0069】

次に、計算粒子のそれぞれについて式（１４）～（１９）の支配方程式を数値的に解くことにより、計算粒子のそれぞれに付与される物理量を計算し、計算粒子のそれぞれの物理量及び位置を更新する（ステップＳ５）。その後、計算粒子のそれぞれについて、その計算粒子の位置における流路の特徴的な寸法に対する空間解像度（計算粒子の直径）を評価する（ステップＳ６）。具体的には、式（２３）、式（２４）が成り立つか否かを判定する。式（２３）が成り立つ場合は空間解像度が不足していると判定し、式（２４）が成り立つ場合は空間解像度が過剰であると判定する。

【0070】

流路の特徴的な寸法に対して空間解像度が不足している場合は、計算粒子を分割する（ステップＳ７）。計算粒子の分割は、例えば図３Ａ～図４Ｂを参照して説明した方法で行う。なお、分割数は、式（２３）のカーネル幅ｈ_ｉと流路の特徴液な寸法ｗ_ｉとの比に基づいて決定するとよい。さらに、流路の特徴的な寸法に対して空間解像度が過剰である場合は、計算粒子を合体させる（ステップＳ８）。計算粒子の合体は、例えば図５Ａ～図５Ｂを参照して説明した方法で行う。なお、ステップＳ８をステップＳ７の前に行ってもよい。

【0071】

ステップＳ５からステップＳ８までの手順を、計算終了条件が満たされるまで繰り返す（ステップＳ９）。計算終了条件が満たされる場合は、計算結果を出力部３２に出力する（ステップＳ１０）。

【0072】

次に、図８Ａ～図９を参照して、本実施例によるシミュレーション方法を用いて流路内の流体の流れをシミュレーションした結果について説明する。解析対象の流路は、図１Ａに示したものと同一の二次元の流路とした。流入端における流路幅Ｗｉは４０ｍｍ、流出端における流路幅Ｗｏは１８ｍｍ、スロート部における流路幅Ｗｔは１０ｍｍ、流路の長さＬは２００ｍｍである。

【0073】

流路の断面形状は、C. A. Moses, G. D. Stein, “On the Growth of Steam Droplets Formed in a Laval Nozzle Using Both Static Pressure and Light Scattering Measurement”, Journal of Fluids Engineering, 100(3), 311-322, September 1987, ISSN 0098-2202, doi: 10.1115/1.3448672, https://doi.org/10.1115/1,3448672に記載された実験で用いられたものと同一である。流れの条件には、上記論文に記載された実験にならって、よどみ点圧力ｐ_０を６７．６６ｋＰａとし、よどみ点温度Ｔ_０を３７６．７Ｋとした。

【0074】

シミュレーションの初期条件において、流れ場全体にカーネル幅ｈ＝４／（２^１／２）ｍｍの計算粒子を、４ｍｍ間隔で配置した。また、計算粒子の各々の圧力及び温度として、それぞれよどみ点圧力及びよどみ点温度を与え、初期速度としてゼロを与えた。

【0075】

図８Ａは、計算開始時点における計算粒子２０の分布を示す図である。計算の時間刻み幅を１．０×１０^－６秒とし、初期状態から流れ場が定常状態に達するまで、計算粒子の位置及び速度を時間発展させた。

【0076】

図８Ｂは、定常状態における計算粒子２０の分布を示す図である。図８Ｂにおいて、計算粒子２０の濃淡は、カーネル幅を表す。流路の特徴的な寸法が相対的に大きい領域Ａ１において、カーネル幅が相対的に大きくなっており、流路の特徴的な寸法が相対的に小さい領域Ａ３において、カーネル幅が相対的に小さくなっている。領域Ａ２、Ａ４におけるカーネル幅は、領域Ａ１におけるカーネル幅と領域Ａ３におけるカーネル幅との中間の大きさになっている。

【0077】

図９は、定常状態における流路の中心軸上での圧力及び液滴の平均半径の時間平均値の分布の計算結果を、上記論文に示された実験結果と比較して示すグラフである。横軸は、流路のスロート部からの距離を単位［ｃｍ］で表し、左縦軸はよどみ点圧力ｐ_０で正規化した正規化圧力を表し、右縦軸は液滴の平均半径を単位［ｎｍ］で表す。液滴の平均粒径ｒは、以下の式により計算した。

【数36】

ここで、ρ_ｌは液相の密度である。液相の密度ρ_ｌの値は、シミュレーション条件として与えられる。

【0078】

図９において実線は、本実施例によるシミュレーション方法による計算結果を示し、丸記号は、上記論文に記載された実験結果を示す。計算によって得られた結果は、圧力分布及び液滴の平均半径の分布ともに、実験結果とほぼ一致していることがわかる。図９に示した計算結果から、本実施例によるシミュレーション方法によって、湿り蒸気流れの解析を行うことが可能であることが確認された。

【0079】

次に、本実施例の優れた効果について説明する。
本実施例では、ステップＳ７、Ｓ８（図７）において、流路の特徴的な寸法に応じて必要かつ十分な空間解像度が自動的に決定される。オペレータが、高い空間解像度を要する領域と、低い空間解像度で解析を行ってもよい領域とを定義する必要がないため、オペレータに、流路内の場所に応じて好適な空間解像度を決定する高い技量が要求されない。

【0080】

また、解析を行う前に、必要とされる空間解像度に応じて、解析空間を複数の領域に区分するといった前処理を行う必要がない。このため、解析に必要とされる工程数を削減することができる。

【0081】

次に、図１０Ａ及び図１０Ｂを参照して、上記実施例の変形例について説明する。上記実施例では、二次元の流路内の流れ場を解析対象としたが、図１０Ａおよび図１０Ｂを参照して説明する変形例では、三次元の流路内の流れ場の解析を行う。図１０Ａ及び図１０Ｂは、変形例によるシミュレーション方法で解析対象となる流路の、流れ方向に直交する断面（以下、単に「流路断面」という。）を示す図である。三次元の流路においても、図１Ｂを参照して説明したように符号付距離関数の値を計算し、図１Ｃを参照して説明したように流路の特徴的な寸法を計算すればよい。

【0082】

図１０Ａに示した変形例においては、流路１０が中心軸Ｃの周りに回転対称であり、流れ方向の任意の位置における流路断面が円形である。流れ方向のある位置における流路断面の直径をＤと標記する。流路断面が円形の流路に定義される符号付距離関数の等値線は同心円となる。また、符号付距離関数の値は、半径方向に関して中心軸Ｃにおいて極大値をとる。このため、任意の点Ｐから勾配ベクトルの方向に、符号付距離関数の極大値を探索すると、流路断面の中心軸Ｃの位置における符号付距離関数の値が極大値として検出される。この極大値の２倍が、点Ｐの位置における流路の特徴的な寸法として採用される。

【0083】

図１０Ｂに示した変形例においては、流路１０の流れ方向の任意の位置における流路断面が正三角形であり、流路断面を流れ方向に移動させたとき流路断面の重心は直線に沿って移動する。流路断面の重心が辿る直線を流路の中心軸Ｃということとする。任意の点Ｐから勾配ベクトルの方向に、符号付距離関数の極大値を探索すると、流れ方向の任意の位置における流路断面の各頂点の二等分線で構成される面（以下、「二等分面Ｂ」という。）と交差する箇所Ｐｍにおいて、極大値が現れる。この極大値の２倍が、点Ｐにおける流路の特徴的な寸法として採用される。

【0084】

このように、三次元の流路についても、流路内の任意の点における流路の特徴的な寸法を定義することができる。例えば、流路の特徴的な寸法は、計算粒子の位置における流路の広がりの程度に依存する量であるということができる。

【0085】

例えば、図１０Ｂに示した正三角形の流路断面を有する流路において、１つの辺の中点の近傍においては、最も近い流路壁に向かう向きには、流路の広がりの程度が小さいが、この流路壁から遠ざかる向き、及びこの流路壁に平行な２方向には、流路の広がりの程度が大きい。流路断面の頂点の近傍においては、最も近い流路壁に向かう向き、この流路壁に平行で頂点の側に向かう向き、及びこの流路壁から遠ざかる向きの３方向に、流路の広がりの程度が小さく、この流路壁に平行で頂点とは反対側に向かう向きの１方向にのみ、流路の広がりの程度が大きい。頂点近傍における流路の広がりの程度は、各辺の中点の近傍における流路の広がりの程度より小さいといえる。

【0086】

流路の空間的な広がりの程度が相対的に小さい頂点近傍における流路の特徴的な寸法が、流路の空間的な広がりの程度が相対的に大きい各辺の中点近傍における流路の特徴的な寸法より小さい。このように、任意の位置における流路の特徴的な寸法は、その位置における流路の空間的な広がりの程度を表す指標ということができる。

【0087】

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

【符号の説明】

【0088】

１０ 流路
２０、２０ａ、２０ａ１、２０ａ２、２０ａｂ、２０ｂ、２０ｃ、２０ｄ、２０ｅ、２０ｓ 計算粒子
３０ 入力部
３１ 処理部
３２ 出力部
３３ 記憶部

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】