7648516 シミュレーション装置、シミュレーション方法、及びプログラム

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

(12)【公報種別】特許公報(B2)

(11)【特許番号】

(24)【登録日】2025-03-10

(45)【発行日】2025-03-18

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

(51)【国際特許分類】

   G16Z 99/00 20190101AFI20250311BHJP

   G06F 30/25 20200101ALI20250311BHJP

【ＦＩ】

G16Z99/00

G06F30/25

【請求項の数】 6

(21)【出願番号】P 2021210790

(22)【出願日】2021-12-24

(65)【公開番号】P2023095105

(43)【公開日】2023-07-06

【審査請求日】2024-01-17

(73)【特許権者】

【識別番号】000002107

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

(74)【代理人】

【識別番号】100105887

【弁理士】

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

(72)【発明者】

【氏名】廣田 真人

【審査官】宮地 匡人

(56)【参考文献】

【文献】特開２００７－１４０９９３（ＪＰ，Ａ）

【文献】特開２０１０－０４９５６１（ＪＰ，Ａ）

【文献】中国特許出願公開第１１３０１１０７５（ＣＮ，Ａ）

【文献】伊澤 精一郎，解像度の異なる小領域を含んだ粒子法の開発，第32回数値流体力学シンポジウム講演論文集(CD-ROM)，2018年，pp.1-2

【文献】田中 正幸，解像度可変型MPS法，日本計算工学会論文集，2009年01月28日，Vol.2009，paper. 20090001

(58)【調査した分野】(Int.Cl.，ＤＢ名)

Ｇ１６Ｚ ９９／００

Ｇ０６Ｆ ３０／２５

ＪＳＴＰｌｕｓ（ＪＤｒｅａｍＩＩＩ）

ＪＳＴ７５８０（ＪＤｒｅａｍＩＩＩ）

(57)【特許請求の範囲】

【請求項1】

粒子法を用いて流体の流れを解析するシミュレーション装置であって、
シミュレーション条件が入力される入力部と、
前記入力部に入力されたシミュレーション条件に基づいて、粒子法を用いて複数の粒子のそれぞれの位置を時間発展させる処理部と
を備え、
前記処理部は、
解析空間を複数の分割領域に分割し、
前記複数の分割領域の各々に、分割領域ごとに大きさの異なる複数の粒子を配置し、
前記複数の分割領域の各々について複数の粒子の位置を時間発展させる際に、複数の分割領域の各々について、分割領域内の現時点の粒子の体積の和が一定になるように、分割領域に新たな粒子を追加するか否かを判定し、
新たな粒子を追加すると判定された分割領域に新たな粒子を追加するシミュレーション装置。

【請求項2】

前記処理部は、
前記複数の分割領域に新たな粒子を追加するか否かの判定において、
前記複数の分割領域のそれぞれの境界に沿って定義した複数の格子点のそれぞれについて、格子点を通って分割領域に流入する粒子の単位時間当たりの個数を時間発展させ、流入する粒子の個数が判定閾値を超えた時点で、流入する粒子の個数が前記判定閾値を超えた格子点に粒子を追加する請求項１に記載のシミュレーション装置。

【請求項3】

前記処理部は、
分割領域に新たな粒子を追加するか否かの判定において、前記複数の格子点の各々について、粒子が分割領域に流入しているか流出しているかを判定し、流出入の判定結果を加味して、前記複数の格子点のそれぞれに新たな粒子を追加するか否かを判定する請求項２に記載のシミュレーション装置。

【請求項4】

前記処理部は、
前記複数の分割領域の各々について粒子法を用いて粒子の挙動を解析する際に、前記複数の分割領域のうち一つの分割領域に着目したとき、着目する分割領域を隣接する分割領域内に向かって拡張したバッファ領域を定義し、着目する分割領域に配置された粒子と同じ大きさの粒子を前記バッファ領域に配置し、前記バッファ領域内の粒子に付与される物理量を、前記バッファ領域が含まれる隣接する分割領域内の粒子に付与された物理量に基づいて計算する請求項１乃至３のいずれか１項に記載のシミュレーション装置。

【請求項5】

粒子法を用いて流体の流れを解析するシミュレーション方法であって、
解析空間を複数の分割領域に分割し、
前記複数の分割領域の各々に、分割領域ごとに大きさの異なる複数の粒子を配置し、
前記複数の分割領域の各々について複数の粒子の位置を時間発展させる際に、複数の分割領域の各々について、分割領域内の現時点の粒子の体積の和が一定になるように、分割領域に新たな粒子を追加するか否かを判定し、
新たな粒子を追加すると判定された分割領域に新たな粒子を追加するシミュレーション方法。

【請求項6】

解析空間を複数の分割領域に分割する手順と、
前記複数の分割領域の各々に、分割領域ごとに大きさの異なる複数の粒子を配置する手順と、
前記複数の分割領域の各々について複数の粒子の位置を時間発展させる際に、複数の分割領域の各々について、分割領域内の現時点の粒子の体積の和が一定になるように、分割領域に新たな粒子を追加するか否かを判定する手順と、
新たな粒子を追加すると判定された分割領域に新たな粒子を追加する手順と
をコンピュータに実行させるプログラム。

【発明の詳細な説明】

【技術分野】

【0001】

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

【背景技術】

【0002】

流体の運動を粒子系の運動として近似して、流体の挙動を解析するシミュレーション方法が知られている。このシミュレーション方法は、粒子法と呼ばれる。粒子法では、流体が複数の粒子で表され、流体の流れ場の中に配置される壁境界も一般的に複数の粒子で表される。

【先行技術文献】

【非特許文献】

【0003】

【文献】原田隆宏、越塚誠一、「ＳＰＨにおける壁境界計算手法の改良」、情報処理学会論文誌、Ｖｏｌ．４８、Ｎｏ．４、１８３８～１８４６頁、２００７年

【発明の概要】

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

【0004】

粒子法による流体解析においては、空間解像度が粒子の直径によって決まる。粒子の直径は、解析しようとする流体運動の長さスケールに対して十分小さく設定しなければならない。一般に、流体運動の長さスケールは場所によって異なる。例えば、流体内に配置される物体の壁から遠い領域では、粒子の直径を大きくしても十分精度の高い解析を行うことができる。ところが、壁の近傍では、流体運動の長さスケールは境界層の厚さで代表されるため、壁近傍の流体の運動を精度よく解析するためには、粒子の直径を境界層の厚さよりも小さくしなければならない。

【0005】

解析空間全体において高精度の解析を行うためには、壁から遠い領域においても、境界層の厚さに応じて粒子の直径を小さくしなければならない。このため、解析空間に配置する粒子数が多くなり、その結果、計算量が増大してしまう。

【0006】

本発明の目的は、粒子法を用いた解析において、解析精度を高く維持し、かつ計算量を削減することが可能なシミュレーション装置、シミュレーション方法、及びプログラムを提供することである。

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

【0007】

本発明の一観点によると、
粒子法を用いて流体の流れを解析するシミュレーション装置であって、
シミュレーション条件が入力される入力部と、
前記入力部に入力されたシミュレーション条件に基づいて、粒子法を用いて複数の粒子のそれぞれの位置を時間発展させる処理部と
を備え、
前記処理部は、
解析空間を複数の分割領域に分割し、
前記複数の分割領域の各々に、分割領域ごとに大きさの異なる複数の粒子を配置し、
前記複数の分割領域の各々について複数の粒子の位置を時間発展させる際に、複数の分割領域の各々について、分割領域内の現時点の粒子の体積の和が一定になるように、分割領域に新たな粒子を追加するか否かを判定し、
新たな粒子を追加すると判定された分割領域に新たな粒子を追加するシミュレーション装置が提供される。

【0008】

本発明の他の観点によると、
粒子法を用いて流体の流れを解析するシミュレーション方法であって、
解析空間を複数の分割領域に分割し、
前記複数の分割領域の各々に、分割領域ごとに大きさの異なる複数の粒子を配置し、
前記複数の分割領域の各々について複数の粒子の位置を時間発展させる際に、複数の分割領域の各々について、分割領域内の現時点の粒子の体積の和が一定になるように、分割領域に新たな粒子を追加するか否かを判定し、
新たな粒子を追加すると判定された分割領域に新たな粒子を追加するシミュレーション方法が提供される。

【0009】

本発明のさらに他の観点によると、
解析空間を複数の分割領域に分割する手順と、
前記複数の分割領域の各々に、分割領域ごとに大きさの異なる複数の粒子を配置する手順と、
前記複数の分割領域の各々について複数の粒子の位置を時間発展させる際に、複数の分割領域の各々について、分割領域内の現時点の粒子の体積の和が一定になるように、分割領域に新たな粒子を追加するか否かを判定する手順と、
新たな粒子を追加すると判定された分割領域に新たな粒子を追加する手順と
をコンピュータに実行させるプログラムが提供される。

【発明の効果】

【0010】

分割領域ごとに大きさの異なる複数の粒子を配置することにより、分割領域ごとに空間解像度を異ならせることができる。高い空間解像度が要求される分割領域の粒子を小さくすることにより、要求される高い空間解像度で解析を行うことができる。高い空間解像度が要求されない分割領域の粒子を大きくすることにより、粒子数を削減することができる。これにより、シミュレーションにおける計算量を削減することができる。分割領域内の現時点の粒子の体積の和が一定になるように新たな粒子を追加することにより、分割領域内の粒子数が変動するような圧縮性の流れ場をシミュレーションすることができる。

【図面の簡単な説明】

【0011】

【図1】図１は、本実施例によるシミュレーション装置の解析対象となる解析モデルの一例を示す斜視図である。

【図2】図２は、解析空間の一部の領域を二次元で表した模式図である。

【図3】図３Ａは、相互に隣接した２つの分割領域の一部分を示す模式図であり、図３Ｂは、２つの分割領域の間で物理量の交換を行う様子を示す模式図である。

【図4】図４は、解析空間の一部を二次元的に示す模式図である。

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

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

【図7】図７Ａ及び図７Ｂは、解析空間に配置した粒子の初期状態を示す図である。

【図8】図８Ａ及び図８Ｂに、流れ場が定常状態に達した状態における粒子分布及びマッハ数分布の瞬時値を示す。

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

【0012】

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

【0013】

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

【数1】

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

【数2】

【0014】

ここで、ｘ_ｉｊは位置ベクトルｘ_ｉと位置ベクトルｘ_ｊとの差の絶対値であり、ｑ＝ｘ_ｉｊ／ｈである。ｈは、平滑化長さであり、カーネル幅ともいわれる。カーネル幅ｈは、複数の粒子を同一の質量及び物性密度の球体に置き換えた等価球の直径に等しい。本明細書において、カーネル幅ｈを粒子の直径という場合がある。

【0015】

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

【数3】

ここで、Ｖ_ｊ＝ｍ_ｊ／ρ_ｊである。下付き文字ｊは、ｊ番目の粒子の持つ物理量の値であることを意味する。

【0016】

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

【数4】

【数5】

【0017】

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

【数6】

【0018】

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

【数7】

【0019】

次に、圧縮性流れを解析対象とする場合の支配方程式について説明する。非粘性の圧縮性流れの支配方程式をＳＰＨ法によって離散化した例を以下に示す。

【数8】

【数9】

【数10】

ここで、ρ_ｉ、ｖ_ｉ、ｐ_ｉ、ｕ_ｉは、それぞれはｉ番目の粒子の密度、速度ベクトル、圧力、及び内部エネルギであり、ｔは時間である。式（８）は質量保存則、式（９）は運動量保存則、式（１０）はエネルギ保存則を表す。

【0020】

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

【0021】

ＳＰＨ法をはじめとする粒子法を用いた流体解析においては、解析の空間解像度は計算粒子の直径によって決まる。計算粒子の直径は、解析しようとする流体運動の長さスケールに対して十分小さく設定しなければならない。一般に、流体運動の長さスケールは場所によって異なる。例えば、物体の壁近傍では流体運動の長さスケールは境界層厚さで代表され、高い空間解像度が必要となる。一方で、壁から離れた領域では、流体運動の長さスケールは境界層厚さよりも大きいことが多い。そのため、壁から離れた領域では、壁近傍の領域に比べて低い空間解像度で解析を行うことができる。数値解析は限られた計算資源で行う必要があるため、高い空間解像度を要する領域では計算粒子径を小さく設定して必要な空間解像度を確保すると同時に、そうでない領域では計算粒子径を大きく設定して解析に要する計算粒子数を節約することが必要である。

【0022】

従来、解析空間を複数の直方体領域に分割し、分割された領域ごとに空間解像度を設定することを可能とする手法が提案されている。しかしながら、壁の形状が複雑な流れ場を解析する場合、直方体のような単純な形状で領域を分割しても、高い空間解像度が必要な領域とそうでない領域とを効果的に分割することは困難である。また、従来手法は非圧縮性流れを対象としており、気体の急激な膨張や収縮、衝撃波の発生を伴う圧縮性の強い流れにおいては用いることができない。

【0023】

図１は、本実施例によるシミュレーション装置の解析対象となる解析モデルの一例を示す斜視図である。解析モデルは、解析空間１０と、解析空間１０に配置された壁１１とを含む。解析空間１０が、複数の領域（以下、分割領域Ωという。）に分割されている。図１に示した例では、解析空間１０が３つの分割領域Ω_１、Ω_２、Ω_３に分割されている。最も内側の第１の分割領域Ω_１が、壁１１を包含しており、第２の分割領域Ω_２が第１の分割領域を包含している。解析空間１０のうち第２の分割領域Ω_２の外側の領域が第３の分割領域Ω_３に相当する。

【0024】

分割領域Ωのそれぞれに複数の粒子２０が配置される。なお、図１では、分割領域Ωのそれぞれに対して１つの粒子２０のみを示している。粒子２０の直径は、分割領域Ωごとに異なっている。例えば、最も内側の第１の分割領域Ω_１内の粒子２０が最も小さく、最も外側の第３の分割領域Ω_３内の粒子２０が最も大きい。１つの分割領域Ω内の複数の粒子２０の大きさは同一である。粒子２０の大きさは、各分割領域Ω内の流体の運動を解析するために必要な空間解像度に応じて決定される。

【0025】

ＳＰＨ法を適用するに際し、解析空間１０を複数の任意形状の分割領域Ωに分割し、分割領域Ωごとに空間解像度を変えるためには、以下の３つの技術要素が必要である。すなわち、第１に、解析空間１０を複数の任意形状の分割領域Ωに分割する技術、第２に、分割領域Ωのうち隣接する領域間で物理量を交換する技術、第３に、分割領域Ωのそれぞれに流入する粒子及びそれぞれから流出する粒子を制御する技術である。以下、これらの技術について説明する。

【0026】

［解析空間を複数の分割領域に分割する技術］
まず、解析空間を複数の分割領域Ωに分割する技術について説明する。
解析空間１０（図１）を任意形状の複数の分割領域Ωに分割する手法に、例えば符号付き距離関数（ＳＤＦ）を用いる。符号付き距離関数φは、解析空間１０の全体において定義される連続関数であり、与えられた任意の境界Ｓから任意の位置ｘまでの最短距離を関数の値として持つ。また、符号付き距離関数φの符号は、境界Ｓの内側では負であり、外側では正である。

【0027】

符号付き距離関数φの値を計算するためには、初期条件として境界Ｓ上でφ＝０を与え、下記の移流方程式を解けばよい。

【数11】

ここで、ｔは時間、ベクトルｖは人工的な移流速度である。

【0028】

解析空間１０を分割するための境界Ｓの位置及び形状を決定し、その境界Ｓを用いて符号付き距離関数φを計算する。符号付き距離関数φの値が等しくなる少なくとも一つの等数値面によって、解析空間１０を複数の分割領域Ωに分割する。等数値面が、相互に接する２つの分割領域Ωの境界面となる。境界Ｓの形状として、例えば解析空間１０内に配置される壁１１の形状に基づいた形状を採用するとよい。

【0029】

境界Ｓの形状は、ＣＡＤデータ等を用いてユーザが任意に決めることができる。このため、解析空間１０を任意の形状の複数の分割領域Ωに分割することができる。

【0030】

次に、図２を参照して、解析空間１０内の離散化された複数の格子点において、符号付き距離関数φの値を求める具体例について説明する。

【0031】

図２は、解析空間１０の一部の領域を二次元で表した模式図である。解析空間１０内に境界Ｓを設定する。さらに、解析空間１０内に直交等間隔格子を張る。境界Ｓの近傍の複数の格子点Ｐのそれぞれについて、格子点Ｐから境界Ｓまでの距離Ｌｐｓの厳密な値を計算する。この値が、格子点Ｐの位置における符号付き距離関数φの値に相当する。境界Ｓの内側の格子点Ｐについては、距離Ｌｐｓとして負の値を与える。境界Ｓの近傍の格子点Ｐ以外の格子点Ｐには、境界Ｓの近傍の格子点Ｐから境界Ｓまでの距離Ｌｐｓより大きな適当な値を初期値として設定する。このとき、境界Ｓの内側の格子点Ｐには負の値を設定し、外側の格子点Ｐには正の値を設定する。

【0032】

上述の式（１１）を離散化して符号付き距離関数φの値を時間発展させる。全ての格子点Ｐにおける符号付き距離関数φの値が収束したら計算を終了する。この収束値が、各格子点Ｐにおける符号付き距離関数φの値に相当する。このように、解析空間１０内に定義した複数の格子点Ｐのそれぞれの位置における符号付き距離関数φの値を決定することができる。符号付き距離関数φの値がほぼ同一の複数の格子点Ｐを含む面を、分割領域Ωの境界として採用する。

【0033】

［隣接する２つの分割領域の間で物理量を交換する技術］
次に、図３Ａ及び図３Ｂを参照して、隣接する２つの分割領域Ωの間で物理量を交換する技術について説明する。

【0034】

図３Ａは、相互に隣接した２つの分割領域Ω_１及び分割領域Ω_２の一部分を示す模式図である。一方の分割領域Ω_１内に複数の粒子２０Ａが配置されており、他方の分割領域Ω_２内に複数の粒子２０Ｂが配置されている。一方の分割領域Ω_１内の粒子２０Ａは、他方の分割領域Ω_２内の粒子２０Ｂより大きい。分割領域Ω_１内の粒子２０Ａの直径をｈ_１と標記し、分割領域Ω_２内の粒子２０Ｂの直径をｈ_２と標記する。分割領域Ω_１内の粒子２０Ａに付す通し番号をｉと標記し、分割領域Ω_２内の粒子２０Ｂに付す通し番号をｊと標記する。

【0035】

図３Ｂは、２つの分割領域Ω_１、Ω_２の間で物理量の交換を行う様子を示す模式図である。２つの分割領域Ω_１、Ω_２の境界面に、バッファ領域Ω_ｂ１、Ω_ｂ２を定義する。バッファ領域Ω_ｂ１は、一方の分割領域Ω_１を他方の分割領域Ω_２内に向かって拡張した領域であり、バッファ領域Ω_ｂ２は、一方の分割領域Ω_２を他方の分割領域Ω_１内に向かって拡張した領域である。分割領域Ω_１を拡張したバッファ領域Ω_ｂ１は、分割領域Ω_２の一部分と重なる。同様に、分割領域Ω_２を拡張したバッファ領域Ω_ｂ２は、分割領域Ω_１の一部分と重なる。

【0036】

バッファ領域Ω_ｂ１、Ω_ｂ２の幅Δｘは、例えば分割領域Ω_１内の粒子２０Ａの直径ｈ_１及び分割領域Ω_２内の粒子２０Ｂの直径ｈ_２のうち大きい方の直径の４倍程度にするとよい。バッファ領域Ω_ｂ１内の粒子２０Ａの物理量Ａ_ｉは、支配方程式を解いて算出する代わりに、バッファ領域Ω_ｂ１に重なっている分割領域Ω_２内の粒子２０Ｂが持つ物理量のカーネル補間によって算出する。すなわち、物理量Ａ_ｉは、以下の式を用いて算出する。

【数12】

式（１２）の右辺のＡ_ｊは、ＳＰＨ法の支配方程式を解いて算出される物理量であり、Σは、分割領域Ω_２内の粒子についての和を意味する。同様に、バッファ領域Ω_ｂ２内の粒子２０Ｂの物理量Ａ_ｊは、以下の式を用いて算出する。

【数13】

式（１３）の右辺のＡ_ｉは、ＳＰＨ法の支配方程式を解いて算出される物理量であり、Σは、分割領域Ω_１内の粒子についての和を意味する。

【0037】

［分割領域に流入する粒子及び分割領域から流出する粒子を制御する技術］
次に、図４及び図５を参照して、分割領域Ωに流入する粒子及び分割領域Ωから流出する粒子を制御する技術について説明する。

【0038】

ＳＰＨ法によるシミュレーションにおいては、流れとともに複数の粒子が移動するため、分割領域Ωから流出する粒子及び分割領域Ωに流入する粒子に合わせて、粒子を除去したり追加したりする処理が必要である。粒子の除去及び追加の処理と流れ場の物理量との整合性が取れていないと、計算の破綻や解析精度の低下が生じる。

【0039】

解析空間１０（図１）を任意の形状に分割した場合、分割領域Ωの境界面における流れの向き及び速さを事前に決めることができない。さらに、流れが非定常である場合には、分割領域Ωの境界面における流れの向き及び速さが時々刻々と変化する。したがって、任意の形状の複数の分割領域Ωに流入する粒子、及び分割領域Ωから流出する粒子を、流れ場の物理量と整合性が取れるように取り扱う手法が必要である。

【0040】

また、流体の持つ圧縮性が強く現れる流れ場においては、流れ場の局所的な密度が時々刻々と変動するため、非圧縮性の流れ場を前提とした制御技術を適用することができない。

【0041】

図４は、解析空間１０の一部を二次元的に示す模式図である。解析空間１０内に２つの分割領域Ω_１、Ω_２が定義されている。分割領域Ω_１から隣接する分割領域Ω_２内に向かってバッファ領域Ω_ｂ１が広がっている。バッファ領域Ω_ｂ１の外側の表面をＳ_１と標記する。

【0042】

バッファ領域Ω_ｂ１の表面Ｓ_１の内側に沿って複数の格子点Ｐ_ｋを配置する。これらの格子点Ｐ_ｋは、例えば解析空間１０の全体に直交等間隔格子を張り、そのうちバッファ領域Ω_ｂ１の内部にあって、表面Ｓ_１からの距離が２^１／２×ｈ_１以内にある格子点を選ぶことで配置できる。ここで、ｈ_１はバッファ領域Ω_１内の粒子の直径（カーネル幅）である。格子点Ｐ_ｋの間隔は、カーネル幅ｈ_１と同じ程度とすることができる。このとき、各格子点上の物理量は、表面Ｓ_１上の物理量と見なすことができる。格子点Ｐ_ｋ上の任意の物理量Ａ_ｋは、カーネル関数の補間公式を用いて以下の式で表すことができる。ここで、ｋは格子点Ｐ_ｋに付された通し番号である。

【数14】

ここで、式（１４）の右辺のΣは、分割領域Ω_２内の粒子について合計することを意味する。

【0043】

式（１４）は、分割領域Ω_１から広がるバッファ領域Ω_ｂ１の表面に沿う格子点Ｐ_ｋの物理量を、分割領域Ω_１に隣接する分割領域Ω_２内部の物理量を補間して算出することを意味する。ここで、補間される物理量Ａｋは、密度ρ_ｋと速度ベクトルｖ_ｋである。

【0044】

格子点の位置における表面Ｓ_１の法線ベクトルｎ_ｋと速度ベクトルｖ_ｋとを用いて、その内積ｎ_ｋ・ｖ_ｋの符号を調べることで、表面Ｓ_１の格子点Ｐ_ｋにおいて流体が流入している部分と流出している部分とを判別することができる。すなわち、内積が負である格子点Ｐ_ｋの位置においては表面Ｓ１を横切ってバッファ領域Ω_ｂ１内に粒子が流入しており、内積が正である格子点Ｐ_ｋの位置においては、バッファ領域Ω_ｂ１から粒子が流出している。このときの体積流量は、（ｎ_ｋ・ｖ_ｋ）ｈ_１ ^２と見積もることができる。

【0045】

内積が負である格子点Ｐ_ｋにおいて、単位時間当たりにバッファ領域Ω_ｂ１の表面を横切って分割領域Ω_１に流入する粒子の個数ｄＮ_ｋ／ｄｔは、以下の式で表すことができる。

【数15】

ここで、ｍ_１は分割領域Ω_１の粒子１個あたりの質量である。式（１５）の分母は、格子点Ｐ_ｋの位置に存在する粒子１個の体積に相当する。式（１５）を積分することにより、格子点Ｐｋごとに、その格子点Ｐ_ｋを通過して分割領域Ωに流入した粒子の個数Ｎ_ｋを求めることができる。

【0046】

分割領域Ω_１に流入させる粒子の位置、個数及び時刻は、式（１５）で与えられる微分方程式を解くことで算出することができるが、流れ場の離散化に伴う誤差により、流入させる個数に不整合が生じてしまう。次に、その不整合を補正する手法について説明する。

【0047】

分割領域Ω_１の体積をＶ_１と表記する。分割領域Ω_１はシミュレーションの前処理の段階で定義されるため、体積Ｖ_１はシミュレーションを通して一定である。流れ場の連続性により、解析空間１０は常に流体によって満たされているという前提条件を用いると、分割領域Ω_１に存在する粒子の体積ｍ_ｉ／ρ_ｉの総和は分割領域Ω_１の体積Ｖ_１と常に一致していなければならない。すなわち、任意の時刻ｔにおいて、以下の式が満たされなければならない。

【数16】

ここで、式（１６）の左辺のΣは、分割領域Ω_１内の粒子について合計することを意味する。

【0048】

式（１６）から、時刻ｔにおける検査体積Ｖ_０に対する粒子の体積の過不足Ｖ_ｉｎ（ｔ）は、以下の式で算出することができる。検査体積Ｖ_０として、初期状態における分割領域Ω_１の体積を採用するとよい。

【数17】

【0049】

式（１５）と式（１７）とを組み合わせて、各格子点Ｐ_ｋに対して、流体粒子を流入させる時刻を決める微分方程式は、最終的に次のように表される。

【数18】

ここで、αは粒子の体積の過不足に基づく補正項の大きさを調整する定数であり、例えば１を選ぶことができる。各格子点Ｐ_ｋ上で流入する粒子の個数Ｎ_ｋを定義し、初期値を０として式（１８）にしたがって流入する粒子の個数Ｎ_ｋを時間発展させる。流入する粒子の個数Ｎ_ｋの値が判定閾値を超えた時刻において、格子点Ｐ_ｋの位置に粒子を追加する。粒子を追加した格子点Ｐ_ｋについては、流入する粒子の個数Ｎ_ｋの値を０に初期設定する。判定閾値として、例えば１を用いる。この操作を繰り返すことで、流れ場の変動に応じて自動的に分割領域Ω_１に粒子が追加される。なお、分割領域Ω_１から流出した粒子は、解析から除外する。

【0050】

式（１８）は、格子点Ｐ_ｋを通って粒子が分割領域Ωに流入している場合に、流入する粒子の個数Ｎ_ｋを増大させ、粒子が分割領域Ωから流出している場合には、流入する粒子の個数Ｎ_ｋを増減させないことを意味する。このように、複数の格子点Ｐ_ｋの各々について、粒子が分割領域に流入しているか流出しているかを判定し、流出入の判定結果を加味して、複数の格子点Ｐ_ｋのそれぞれに新たな粒子を追加するか否かを決定する。

【0051】

上述の手法を、解析空間１０内に存在するすべての分割領域Ωに対して適用することによって、空間解像度の異なる流れ場同士を接続することができる。

【0052】

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

【0053】

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

【0054】

［シミュレーション方法］
図６は、本実施例のシミュレーション方法の手順を示すフローチャートである。図６に示した各ステップは、処理部３１（図５）が記憶部３３（図５）に記憶されたプログラムを実行することにより実現される。

【0055】

まず、処理部３１が入力部３０に入力されたシミュレーション条件を取得する（ステップＳ１）。シミュレーション条件には、シミュレーション対象の流体を定義する情報、初期条件、境界条件、解析空間１０を複数の分割領域Ω（図１）に分割するための情報等が含まれる。流体を定義する情報には、流体の密度、粘度等の物性値が含まれる。初期条件には、複数の粒子２０（図１）を解析空間１０に配置するための情報、粒子２０の初期速度を指定する情報等が含まれる。境界条件には、解析空間１０の形状及び大きさを指定する情報が含まれる。解析空間１０を複数の分割領域Ω（図１）に分割するための情報には、分割領域Ωの各々の形状及び大きさを指定する情報、分割領域Ωの各々のカーネル幅を指定する情報が含まれる。

【0056】

処理部３１は、シミュレーション条件を取得すると、解析空間１０を複数の分割領域Ωに分割する（ステップＳ２）。解析空間１０の分割には、例えば図２を参照して説明したように、符号付き距離関数φを用い、符号付き距離関数φの等数値面を、複数の分割領域Ωの境界面とする。

【0057】

処理部３１は、複数の分割領域Ωの各々にバッファ領域を定義する（ステップＳ３）。例えば、図３Ｂに示したように、相互に隣接する２つの分割領域Ω_１、Ω_２の一方の分割領域Ω_１を他方の分割領域Ω_２内に拡張してバッファ領域Ω_ｂ１を定義する。同様に、分割領域Ω_２を分割領域Ω_１内に拡張してバッファ領域Ω_ｂ２を定義する。次に、処理部３１は、分割領域Ω_１、Ω_２及びバッファ領域Ω_ｂ１、Ω_ｂ２等の各々に、分割領域Ωごとに直径（カーネル幅）の異なる複数の粒子２０（図１）を配置する。ここで、粒子２０の大きさは、式（２）のカーネル幅ｈで特定される（ステップＳ４）。

【0058】

処理部３１は、複数の粒子２０を配置した後、バッファ領域Ω_ｂ１、Ω_ｂ２等の粒子２０について、図３Ｂを参照して説明した方法により、バッファ領域Ω_ｂ１、Ω_ｂ２等の粒子２０に付与される物理量を決定する（ステップＳ５）。すべての粒子２０に物理量を付与した後、処理部３１は、分割領域Ωの各々について式（８）～（１０）に示した支配方程式を解くことにより、粒子２０に付与される物理量、すなわち速度及び密度を計算する（ステップＳ６）。

【0059】

次に、処理部３１は、分割領域Ωの各々について粒子数を制御する（ステップＳ７）。粒子数の制御は、図４を参照して説明した手法を適用することにより行う。具体的には、格子点Ｐ_ｋ（図４）ごとに微分方程式（１８）を用いて流入する粒子の個数Ｎ_ｋを時間発展させ、流入する粒子の個数Ｎ_ｋの値が判定閾値に達したら、その格子点Ｐ_ｋの位置に粒子を追加する。追加された粒子に付与する物理量は、格子点Ｐ_ｋが属する分割領域Ω内の粒子の物理量を、カーネル関数を用いて補間することにより決定する。

【0060】

粒子数を制御した後、分割領域Ωの各々の粒子の物理量及び位置を更新する（ステップＳ８）。物理量及び位置を時間発展させる処理を終了するまで、ステップＳ５～Ｓ８を繰り返す。粒子の物理量及び位置を時間発展させる処理を終了させる場合は、解析結果を出力部３２（図５）に出力する（ステップＳ９）。

【0061】

［実施例の優れた効果］
次に、上記実施例の優れた効果について説明する。
上記実施例では、複数の分割領域Ω（図１）ごとに異なる空間解像度、すなわち粒子２０の直径を設定してシミュレーションを行うことができる。このため、高い空間解像度が要求される領域の粒子の直径を小さくして、空間解像度の要求を満たすことができる。また、高い空間解像度が要求されない領域においては、粒子の直径を大きくして、計算量を低減させることができる。このように、必要とされる空間解像度を維持し、かつ計算量を低減させることができる。

【0062】

また、上記実施例では、図２を参照して説明したように、任意の形状の境界Ｓを決定し、符号付き距離関数φを用いて複数の分割領域Ωの境界面を決定することができる。例えば、流体の流れ場の壁面や、流れ場内に配置された物体の表面の近傍で高い空間解像度が要求される。このように、高い空間解像度が要求される領域は、流れ場に配置される壁面の形状に応じて決定することが好ましい。本実施例では、壁面の形状に応じて、図２に示した境界Ｓを設定することにより、壁面の形状が反映された形状の分割領域Ωを設定することができる。例えば、本実施例によるシミュレーション方法は、複雑な形状の物体周りの流れ場や、複雑な形状を持つ流路内の流れ場の解析に適用することが可能である。

【0063】

さらに、上記実施例では、微分方程式（１８）を解いて、分割領域Ωの境界面に沿って配置された格子点Ｐ_ｋごとに分割領域Ωに流入する粒子の個数Ｎ_ｋを求め、分割領域Ωに新たに粒子を追加している。このように、分割領域Ω内の現時点の粒子の体積の和が一定になるように粒子数が制御される。このため、分割領域Ω内の粒子数が変動する圧縮性の強い流れ場の解析を行うことができる。

【0064】

［効果の検証］
上記実施例の優れた効果を検証するために、上記実施例によるシミュレーション方法を用いて、タービンの翼列を模した流れ場のシミュレーションを行った。以下、図７Ａ～図８Ｂを参照してこのシミュレーションについて説明する。

【0065】

図７Ａ及び図７Ｂは、解析空間１０に配置した粒子の初期状態を示す図である。翼弦の長さＬｖを１３５．７ｍｍとし、ピッチＬｙを８７．５９ｍｍとし、翼の軸方向の長さＬｘを１０２．５ｍｍとした。翼の断面形状及び寸法は、「White, A. J., Young, J. B. and Walters, P. T., Experimental validation of condensing flow theory for a stationary cascade of steam turbine blades, Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 354, 1704 (1996), pp. 59-88.」に開示された実験で用いられたものと同一とした。

【0066】

解析空間１０の高さを翼列のピッチＬｙと等しくした。流入境界Ｓｉｎから翼の前縁までの距離、及び翼の後縁から流出境界Ｓｏｕｔまでの距離を、それぞれ翼の軸方向の長さＬｘと等しくした。

【0067】

図７Ａは解析空間１０の全体に対してカーネル幅ｈ_０を１ｍｍに設定した場合を示す。図７Ｂは上記実施例によるシミュレーション方法を適用して、解析空間１０を４個の分割領域Ω_０、Ω_１、Ω_２、Ω_３に分割した場合を示す。分割領域Ω_０、Ω_１、Ω_２は、それぞれ翼面からの距離が４ｍｍ以内の領域、１６ｍｍ以内の領域、４０ｍｍ以内の領域である。分割領域Ω_３は、翼の前縁から上流側に４０ｍｍ離れた位置から、翼の後縁から下流側に４０ｍｍ離れた位置までの領域である。分割領域Ω_０、Ω_１、Ω_２、Ω_３のカーネル幅ｈ_０、ｈ_１、ｈ_２、ｈ_３を、それぞれ１ｍｍ、２ｍｍ、４ｍｍ、８ｍｍに設定した。

【0068】

流入境界Ｓｉｎには一様流入条件を適用し、流出境界Ｓｏｕｔにはノイマン条件を課した。解析空間１０の上下方向には周期境界条件を与えた。壁面には断熱条件及びすべり条件を課した。流れのよどみ点圧力とよどみ点温度とを、それぞれ４１．９ｋＰａ、３７８．５Ｋに設定した。計算の初期条件においては、流れ場全体に、分割領域Ωごとに設定されたカーネル幅の粒子を等間隔に配置し、速度をゼロとして、よどみ点圧力及びよどみ点温度を与えた。解析の時間刻み幅は、共に２．６７×１０^－７秒とした。

【0069】

計算の初期状態における粒子の個数は、図７Ａに示した条件で１９９，５９２個、図７Ｂに示した条件で４７，２１６個であった。すなわち、本実施例では、解析空間１０を分割しない手法と比べて粒子数を約１／４に削減することができた。計算に要する時間も同様に、１／４程度に削減された。

【0070】

図８Ａ及び図８Ｂは、流れ場が定常状態に達した状態における粒子分布及びマッハ数分布の瞬時値を示す図である。マッハ数が相対的に小さい粒子を、相対的に淡い色で表している。図８Ａが解析空間１０を分割しない場合（図７Ａ）、図８Ｂが分割領域Ωに分割した場合（図７Ｂ）の解析結果である。流れは翼列を通過する際に膨張し、超音速になる。図８Ａからわかるように、翼の後縁近傍から、マッハ数が急激に変化する領域が直線状に伸びており、その境界で衝撃波が発生していることがわかる。図８Ｂでは、カーネル幅が比較的小さい分割領域Ω_０、Ω_１に、マッハ数が急激に変化する領域が現れている。

【0071】

このシミュレーション結果から、上記実施例によるシミュレーション方法は、流体の急激な膨張や衝撃波の発生を伴い、流体の圧縮性が強く現れる流れ場のシミュレーションを効率的に行うことが可能であることが確認された。

【0072】

なお、図８Ｂに示したシミュレーション結果では、分割領域Ω_２、Ω_３でマッハ数の変化が緩やかになっているように見える。これは、カーネル幅が大きいことによって数値粘性が増大し、それにより衝撃波がなまされた結果であると考えられる。より詳細なシミュレーションを行う必要がある場合は、分割領域Ωの設定、及びカーネル幅を変化させて再度シミュレーションを行えばよい。

【0073】

上述の実施例は例示であり、本発明は上述の実施例に制限されるものではない。例えば、種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。上記実施例では、粒子法としてＳＰＨ法を用いているが、その他の粒子法、例えば散逸粒子法、ＭＰＳ法等を用いた流体解析にも、上記実施例による手法を適用することが可能である。

【符号の説明】

【0074】

１０ 解析空間
１１ 壁
２０、２０Ａ、２０Ｂ 粒子
３０ 入力部
３１ 処理部
３２ 出力部
３３ 記憶部

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】