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

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

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

(11)【公開番号】P2025004498

(43)【公開日】2025-01-15

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

(51)【国際特許分類】

   G16Z 99/00 20190101AFI20250107BHJP

   G01N 11/00 20060101ALI20250107BHJP

【ＦＩ】

G16Z99/00

G01N11/00 A

【審査請求】未請求

【請求項の数】6

【出願形態】ＯＬ

(21)【出願番号】P 2023104219

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

(71)【出願人】

【識別番号】000002107

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

(74)【代理人】

【識別番号】100105887

【弁理士】

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

(72)【発明者】

【氏名】廣田 真人

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

5L049

【Ｆターム（参考）】

5L049DD02

(57)【要約】

【課題】壁によって相互に隔てられた領域を含む流れ場の解析を、粒子法を用いて高精度に行うことが可能なシミュレーション装置を提供する。
【解決手段】入力部にシミュレーション条件が入力される。処理部が、シミュレーション条件に基づいて流体を複数の計算粒子で表し、複数の計算粒子のそれぞれについて支配方程式を数値的に解いて、複数の計算粒子の位置及び物理量を時間発展させる。支配方程式を解く際に、着目する計算粒子に対して、前記壁面を隔てて配置された他の計算粒子を除外する。壁面から複数の計算粒子のそれぞれの位置までの最短経路の長さ及び向きで定義される最短距離ベクトルを示す情報を算出する。２つの計算粒子のそれぞれの位置ベクトル、及び２つの計算粒子のそれぞれに対する最短距離ベクトルを用いて、２つの計算粒子が相互に壁面を隔てて配置されているか否かを判定する。
【選択図】図４

【特許請求の範囲】

【請求項1】

壁面を有する空間内の流体の流れを、粒子法を用いて解析するシミュレーション装置であって、
解析すべき流体を定義する情報、解析の初期条件と境界条件、解析空間に配置される壁面の形状を定義する情報が入力される入力部と、
前記入力部に入力された情報に基づいて、前記流体を複数の計算粒子で表し、前記複数の計算粒子のそれぞれについて支配方程式を数値的に解いて、前記複数の計算粒子の位置及び物理量を時間発展させる処理部と
を備え、
前記処理部は、
前記支配方程式を解く際に、着目する計算粒子に対して、前記壁面を隔てて配置された他の計算粒子を除外して前記支配方程式を解く機能と、
前記壁面から前記複数の計算粒子のそれぞれの位置までの最短経路の長さ及び向きで定義される最短距離ベクトルを示す情報を算出する機能と、
２つの計算粒子のそれぞれの位置ベクトル、及び２つの計算粒子のそれぞれに対する前記最短距離ベクトルを用いて、２つの計算粒子が相互に前記壁面を隔てて配置されているか否かを判定する機能と
を有するシミュレーション装置。

【請求項2】

前記処理部は、
前記壁面のうち前記最短距離ベクトルを与える領域とは異なる領域から前記複数の計算粒子それぞれの位置までの経路の長さの極小値及び経路の向きで定義される極小距離ベクトルを示す情報を算出する機能を、さらに有し、
２つの計算粒子が相互に前記壁面を隔てて配置されているか否かを判定する際に、２つの計算粒子のそれぞれの位置ベクトル、及び２つの計算粒子のそれぞれに対する前記最短距離ベクトルに加えて、２つの計算粒子のそれぞれに対する前記極小距離ベクトルを用いて、２つの計算粒子が相互に前記壁面を隔てて配置されているか否かを判定する請求項１に記載のシミュレーション装置。

【請求項3】

前記解析空間に張られた格子の複数の格子点のそれぞれの位置において、前記壁面から格子点までの前記最短距離ベクトルを示す情報が格納された記憶部を、さらに備え、
前記処理部は、前記複数の計算粒子のそれぞれの位置において、前記複数の格子点のそれぞれの位置に対応する前記最短距離ベクトルを示す情報を補間することにより、前記複数の計算粒子のそれぞれの位置における前記最短距離ベクトルを求める機能を、さらに有する請求項２に記載のシミュレーション装置。

【請求項4】

前記処理部は、前記複数の格子点のそれぞれの位置において、着目する格子点に対する前記最短距離ベクトル、着目する格子点の位置ベクトル、着目する格子点から所定距離以内に存在する複数の他の格子点である参照格子点のそれぞれの位置ベクトルと着目する格子点の位置ベクトルとの差、及び前記複数の参照格子点のそれぞれに対する前記最短距離ベクトルを用いて、着目する格子点に対する前記極小距離ベクトルを示す情報を算出する機能を、さらに有する請求項２または３に記載のシミュレーション装置。

【請求項5】

壁面を有する空間内の流体の流れを、粒子法を用いて解析するシミュレーション方法であって、
シミュレーション条件を与えられたコンピュータが、与えられたシミュレーション条件に基づいて、
前記流体を複数の計算粒子で表し、前記複数の計算粒子のそれぞれについて支配方程式を数値的に解いて、前記複数の計算粒子の位置及び物理量を時間発展させ、
前記支配方程式を解く際に、着目する計算粒子に対して、前記壁面を隔てて配置された他の計算粒子を除外して前記支配方程式を解き、
前記壁面から前記複数の計算粒子のそれぞれの位置までの最短経路の長さ及び向きで定義される最短距離ベクトルを示す情報を算出し、
２つの計算粒子のそれぞれの位置ベクトル、及び２つの計算粒子のそれぞれに対する前記最短距離ベクトルを用いて、２つの計算粒子が相互に前記壁面を隔てて配置されているか否かを判定するシミュレーション方法。

【請求項6】

壁面を有する空間内の流体の流れを、粒子法を用いて解析する機能をコンピュータに実現させるプログラムであって、
与えられたシミュレーション条件に基づいて、
前記流体を複数の計算粒子で表し、前記複数の計算粒子のそれぞれについて支配方程式を数値的に解いて、前記複数の計算粒子の位置及び物理量を時間発展させる機能と、
前記支配方程式を解く際に、着目する計算粒子に対して、前記壁面を隔てて配置された他の計算粒子を除外して前記支配方程式を解く機能と、
前記壁面から前記複数の計算粒子のそれぞれの位置までの最短経路の長さ及び向きで定義される最短距離ベクトルを示す情報を算出する機能と、
２つの計算粒子のそれぞれの位置ベクトル、及び２つの計算粒子のそれぞれに対する前記最短距離ベクトルを用いて、２つの計算粒子が相互に前記壁面を隔てて配置されているか否かを判定する機能と
を実現させるプログラム。

【発明の詳細な説明】

【技術分野】

【0001】

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

【背景技術】

【0002】

流体解析技術として、連続な場を仮想的な粒子（以下、計算粒子という。）の集合を用いて離散化して支配方程式を解くことにより、流れ場の時間発展を解析する粒子法が知られている。与えられた物体周りの流れを解析する場合には、解析条件に壁面の境界条件が課される。例えば、特許文献１及び特許文献２に、任意形状の物体に対して壁面境界条件を課す手法が記載されている。

【先行技術文献】

【特許文献】

【0003】

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

【特許文献2】特開２０２２－８０６４３号公報

【発明の概要】

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

【0004】

粒子法を用いて支配方程式を解く場合には、着目する計算粒子の近傍にある他の計算粒子を参照し、それらに付与された物理量を用いて、着目する計算粒子に付与された物理量を更新する。それぞれの計算粒子に対して設定される参照範囲の内部で、かつ薄い壁によって隔てられた領域に流れが存在する場合、着目する計算粒子の参照範囲が、壁によって隔てられた領域の計算粒子まで及んでしまう。

【0005】

一般に、壁によって相互に隔てられた領域を流れる流体の間では、相互作用はほとんど発生しない。ところが、従来の粒子法を用いて流れ場の解析を行うと、壁によって相互に隔てられた２つの計算粒子の間で物理量のやりとりが発生してしまう。その結果、解析精度が低下してしまう。

【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】

ここで、ｕは速度ベクトル、ｔは時間、ρは密度、ｐは圧力、μは粘性係数を表す。

【0022】

［ＳＰＨ法を用いた支配方程式の離散化手法］
式（９）、式（１０）に示した流れの支配方程式を、ＳＰＨ法を用いて離散化すると、例えば以下のように表すことができる。

【数11】

【数12】

【数13】

【数14】

ここで、下付き文字ｉ及びｊは、計算粒子に付された通し番号を表す。式（１３）は、流体が弱い圧縮性を持つという仮定により導出される状態方程式であり、ｃは仮想的な音速、ρ_０は基準密度を表す。音速ｃの値は、例えば流れ場の代表速度の１０倍程度に設定することができる。基準密度ρ_０の値は、例えば作動流体の密度に設定することができる。式（１１）～式（１４）に示す支配方程式は、非圧縮性流れのＳＰＨ法による解法として広く用いられる。

【0023】

［第１実施例］
次に、第１実施例によるシミュレーション装置について説明する。第１実施例においては、解析空間に任意形状の物体が配置されており、この物体周りの流れを解析する。物体の形状は、例えばＣＡＤデータで与えられる。物体表面（以下、壁面という場合がある。）における境界条件を与えるために、複数の計算粒子のそれぞれの位置から物体表面までの距離を算出する必要がある。

【0024】

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

【数22】

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

【0025】

［解析空間に薄い壁が存在する場合の課題］
次に、図１を参照して、解析空間に薄い壁が存在する場合の課題について説明する。
図１は、解析空間に配置された物体１０の断面図である。物体１０は、流体が収容されている流路を取り囲む外側の壁１０Ａと、流路内に配置された薄い壁１０Ｂとを含む。流路内に複数の計算粒子２０が配置されている。薄い壁１０Ｂの両側にも複数の計算粒子２０が配置される。

【0026】

支配方程式を解く対象となる計算粒子２０を、着目する計算粒子２０ｉということとする。図１において、着目する計算粒子２０ｉに、右下がりの相対的に濃いハッチングを付している。計算粒子２０ｉについて支配方程式を解くときの参照範囲２３ｉは、例えば着目する計算粒子２０ｉの位置を中心とする球形であり、その半径は、例えば計算粒子２０の直径の３倍～１０倍程度である。図１において、参照範囲２３ｉの外周を破線で示している。

【0027】

着目する計算粒子２０ｉが薄い壁１０Ｂの近傍に位置する場合、着目する計算粒子２０ｉに対して薄い壁１０Ｂによって隔てられた反対側の空間まで参照範囲２３ｉが広がる。図１において、参照範囲２３ｉ内に存在する計算粒子２０のうち、薄い壁１０Ｂを境として、着目する計算粒子２０ｉと同じ側に存在する計算粒子２０ａに右上がりの相対的に淡いハッチングを付し、反対側に存在する計算粒子２０ｂを灰色で表している。

【0028】

薄い壁１０Ｂによって隔てられた両側の空間を流れる流体は、相互に影響を及ぼし合わない。このため、着目する計算粒子２０ｉについて支配方程式を解くとき、計算粒子２０ｂを除外することが好ましい。ところが、従来のＳＰＨ法を用いた解析では、着目する計算粒子２０ｉについて支配方程式を解くとき、計算粒子２０ａ及び計算粒子２０ｂの両方が参照されてしまう。以下に説明する第１実施例では、計算粒子２０ｂを除外して、着目する計算粒子２０ｉについて支配方程式を解く。

【0029】

次に、図２Ａ及び図２Ｂを参照して、２つの計算粒子が、相互に薄い壁１０Ｂを隔てて配置されているか否かを判定する方法について説明する。図２Ａは、着目する計算粒子２０ｉと、参照範囲２３ｉ（図１）内の他の計算粒子２０ｊ（参照格子点）とが、薄い壁１０Ｂを隔てて配置されている場合の模式図であり、図２Ｂは、着目する計算粒子２０ｉと、参照範囲２３ｉ（図１）内の他の計算粒子２０ｊとが、薄い壁１０Ｂに関して同じ側に存在する場合の模式図である。着目する計算粒子２０ｉの位置ベクトルをｘ_ｉと標記し、他の計算粒子２０ｊの位置ベクトルをｘ_ｊと標記する。

【0030】

図２Ａにおいては、薄い壁１０Ｂの一方の壁面１０Ｆが、着目する計算粒子２０ｉの側を向き、他方の壁面１０Ｆが、計算粒子２０ｊの側を向く。図２Ｂにおいては、薄い壁１０Ｂの一方の壁面１０Ｆが、着目する計算粒子２０ｉ及び他の計算粒子２０ｊの側を向く。

【0031】

解析空間に直交等間隔格子の複数の格子点が定義されており、複数の格子点のそれぞれの位置において、与えられた物体１０（図１）に対する符号付距離関数が計算され、その値が格子点に関連付けられて記憶されている。

【0032】

まず、着目する計算粒子２０ｉ及び他の計算粒子２０ｊのそれぞれについて、壁面１０Ｆからの最短距離ベクトルｄ_ｉ、ｄ_ｊを算出する。最短距離ベクトルは、壁面１０Ｆから計算粒子２０ｉ、２０ｊのそれぞれの位置までの最短経路の長さ及び向きを持つベクトルである。例えば、計算粒子２０ｉの近傍の複数の格子点のそれぞれに関連付けられて記憶されている符号付距離関数の値、及びその勾配ベクトルに基づいて、線形補間演算を行うことにより、最短距離ベクトルｄ_ｉを算出することができる。最短距離ベクトルｄ_ｊについても同様に算出することができる。

【0033】

計算粒子２０ｉ、２０ｊに対する壁面１０Ｆの位置ベクトルｘ_ｗｉ、ｘ_ｗｊを、以下の式により計算する。

【数23】

【数24】

【0034】

以下の２つの式のいずれかが満たされる場合、計算粒子２０ｉと計算粒子２０ｊとは、相互に壁面１０Ｆを隔てて配置されていると判定することができる。

【数25】

【数26】

ここで、γは任意のパラメータであり、例えばγの値としてゼロを用いることができる。

【0035】

次に、式（２５）及び式（２６）の物理的な意味について説明する。
式（２５）の左辺は、ベクトルｘ_ｉ－ｘ_ｗｊと最短距離ベクトルｄ_ｊとのなす角度θ_ｊｉを表し、式（２６）の左辺は、ベクトルｘ_ｊ－ｘ_ｗｉと最短距離ベクトルｄ_ｉとのなす角度θ_ｉｊを表す。γにゼロを設定した場合、角度θ_ｊｉ、θ_ｉｊの少なくとも一方が９０°より大きい場合、式（２５）及び式（２６）の少なくとも一方が満たされる。

【0036】

図２Ａに示すように、角度θ_ｊｉ、θ_ｉｊの両方が９０°より大きい場合は、計算粒子２０ｉと計算粒子２０ｊとは、相互に壁面１０Ｆを隔てて配置されていると判定される。図２Ｂに示すように、角度θ_ｊｉ、θ_ｉｊの両方が９０°より小さい場合は、計算粒子２０ｉと計算粒子２０ｊとは、相互に壁面１０Ｆによって隔てられていないと判定される。

【0037】

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

【0038】

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

【0039】

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

【0040】

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

【0041】

処理部３１は、与えられた物体形状について、直交等間隔格子の格子点のそれぞれの位置における符号付距離関数の値を計算する（ステップＳ２）。

【0042】

次に、処理部３１は、入力されたシミュレーション条件及び符号付距離関数の値に基づいて、解析空間に複数の計算粒子２０を配置する（ステップＳ３）。具体的には、符号付距離関数の値が正の領域に、複数の計算粒子２０を配置する。その後、複数の計算粒子２０のそれぞれ（着目する計算粒子２０ｉ）について、参照範囲２３ｉ内の計算粒子２０ｊをリスト化した近接粒子リストを作成する（ステップＳ４）。

【0043】

次に、近接粒子リストを作成するときに着目した計算粒子２０ｉと、近接粒子リストに含まれる他の計算粒子２０ｊのそれぞれとが、相互に壁面を隔てて配置されているか否かを判定し、着目する計算粒子２０ｉに対して壁面を隔てて配置されている他の計算粒子２０ｊを近接粒子リストから除外する（ステップＳ５）。着目する計算粒子２０ｉと、近接粒子リストに含まれる他の計算粒子２０ｊのそれぞれとが、相互に壁面を隔てて配置されているか否かの判定には、式（２５）及び式（２６）を用いる。

【0044】

複数の計算粒子のそれぞれについて、式（１１）～式（１４）に示した支配方程式を数値的に解き、計算粒子２０のそれぞれの物理量及び位置を更新する（ステップＳ６）。着目する計算粒子２０ｉについて支配方程式を解くとき、ステップＳ５で一部の計算粒子２０ｊを除外する処理を行った後の近接粒子リストに含まれる複数の計算粒子２０ｊを参照する。すなわち、着目する計算粒子２０に対して、壁面を隔てて配置された他の計算粒子を除外して支配方程式を解く。

【0045】

次に、近接粒子リストを更新するか否かを判定する（ステップＳ７）。例えば、複数の計算粒子２０のそれぞれについて、現時点の近接粒子リストを作成した時点の位置からの移動量が判定上限値を超えたか否かを判定する。少なくとも１つの計算粒子２０の移動量が判定上限値を超えている場合、近接粒子リストを更新すると判定し、それ以外の場合は、近接粒子リストを更新しないと判定する。

【0046】

近接粒子リストを更新すると判定した場合、近接粒子リストを更新する（ステップＳ８）。近接粒子リストの更新は、ステップＳ４で行った近接粒子リストの作成と同一の手順で行う。さらに、更新後の近接粒子リストについて、着目する計算粒子２０ｉに対して壁面を隔てて配置されている他の計算粒子を近接粒子リストから除外する（ステップＳ９）。この手順は、ステップＳ５の手順と同一である。

【0047】

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

【0048】

次に、第１実施例のすぐれた効果について説明する。
第１実施例では、解析空間に薄い壁１０Ｂ（図１）が配置され、その両側を流体が流れる場合であっても、ステップＳ５及びステップＳ９（図４）において、着目する計算粒子２０ｉに対して壁面を隔てて配置されている他の計算粒子を近接粒子リストから除外するため、解析精度の低下を招くことなく解析することができる。

【0049】

さらに、第１実施例では、図２Ａ及び図２Ｂを参照して説明したように、符号付距離関数の値のみに基づいて、２つの計算粒子２０ｉ、２０ｊが、相互に壁面を隔てて配置されているか否かを判定する。このため、解析空間に配置されている壁面の形状及び位置を定義する情報を参照して、壁面を探索する処理を行う必要が無い。

【0050】

次に、図５を参照して第１実施例の変形例について説明する。
図５は、第１実施例の変形例によるシミュレーション方法を説明するための計算粒子２０ｉ、２０ｊ、及び壁面１０Ｆの位置関係を示す模式図である。図２Ａ及び図２Ｂに示した例では、計算粒子２０ｉ、２０ｊのそれぞれの最短距離ベクトルｄ_ｉ、ｄ_ｊが、壁面１０Ｆに対して垂直である。これに対して図５に示した例では、計算粒子２０ｉ、２０ｊが壁面１０Ｆの縁の近傍に配置されている。計算粒子２０ｉの最短距離ベクトルｄ_ｉは、壁面１０Ｆに対して垂直であるが、計算粒子２０ｊの最短距離ベクトルｄ_ｊの始点が壁面１０Ｆの縁に位置し、計算粒子２０ｊの最短距離ベクトルｄ_ｊが壁面１０Ｆに対して垂直ではない。

【0051】

図５に示した例においても、図２Ａに示した例と同様に、角度θ_ｉｊ及び角度θ_ｊｉの両方が９０°より大きい。このため、式（２５）及び式（２６）のパラメータγとしてゼロを用いると、着目する計算粒子２０ｉについて支配方程式を解く際に、計算粒子２０ｊが計算対象から除外される。

【0052】

ところが、図５に示した例では、計算粒子２０ｉの位置における流体の流れと、計算粒子２０ｊの位置における流体の流れとが、相互に影響を及ぼし合うと考えられる。このため、着目する計算粒子２０ｉについて支配方程式を解く際に、計算粒子２０ｊを参照することが好ましい。

【0053】

式（２５）及び式（２６）のパラメータγを、－１より大きく０以下の値に設定すると、角度θ_ｉｊ及びθ_ｊｉが９０°より大きい場合も、計算粒子２０ｉと計算粒子２０ｊとが、相互に壁面１０Ｆによって隔てられていないと判定される。これにより、着目する計算粒子２０ｉについて支配方程式を解く際に、計算粒子２０ｊが参照されることになる。パラメータγの値を調整すると、２つの計算粒子２０ｉ、２０ｊが、相互に壁面１０Ｆを隔てて配置されていると判定する角度θ_ｉｊ、θ_ｊｉの条件が変わる。

【0054】

パラメータγの値は、解析空間に配置されている壁面１０Ｆの形状、流体の特性等に基づいて、適切に設定することが好ましい。

【0055】

また、第１実施例では、ステップＳ４（図４）で近接粒子リストを作成し、ステップＳ８（図４）で近接粒子リストを更新する。近接粒子リストは、支配方程式を解くタイムステップごとに行うのではなく、近接粒子リストを更新する条件を満たす場合だけ、近接粒子リストの更新を行う。この近接粒子リストの更新の際に、計算粒子のそれぞれについて壁面を隔てて配置された他の計算粒子を近接粒子リストから除外する演算を行う。このため、タイムステップごとに、計算粒子のそれぞれについて壁面を隔てて配置された他の計算粒子を近接粒子リストから除外する演算を行う場合と比べて、計算負荷の増大を抑制することができる。

【0056】

［第２実施例］
次に、図６～図１０Ｄを参照して、第２実施例によるシミュレーション装置及びシミュレーション方法について説明する。以下、図１～図４を参照して説明した第１実施例によるシミュレーション装置及びシミュレーション方法と共通の構成については説明を省略する。

【0057】

図６は、第２実施例によるシミュレーション装置の解析対象となる壁面１０Ｆの形状の一例を示す模式図である。解析空間に配置された物体１０の壁面１０Ｆが、一例として、図６に示したように角部分を含むような条件下で着目する計算粒子２０ｉについて支配方程式を解く際に、参照から除外すべき計算粒子２０ｊが除外されない場合があることが判明した。

【0058】

図６に示すように、物体１０が、流路を画定する外側の壁１０Ａと、外側の壁１０Ａから流路内に突き出た薄い壁１０Ｂとを含む。薄い壁１０Ｂの両側を流体が流れる。外側の壁１０Ａと薄い壁１０Ｂとの接続箇所が、壁面１０Ｆの角部を構成する。

【0059】

計算粒子２０ｉ及び計算粒子２０ｊが、相互に薄い壁１０Ｂを隔てて配置されている。ただし、計算粒子２０ｉ、２０ｊに対応する壁面の位置ベクトルｘ_ｗｉ、ｘ_ｗｊで定義される位置は、薄い壁１０Ｂの壁面１０Ｆではなく、外側の壁１０Ａの壁面１０Ｆ上に位置する。２つの計算粒子２０ｉ、２０ｊは、外側の壁１０Ａの壁面１０Ｆから見て同じ側に配置されている。このため、角度θ_ｉｊ、θ_ｊｉは、ともに９０°未満である。図２Ａ及び図２Ｂを参照して説明した方法では、計算粒子２０ｉと計算粒子２０ｊとは、相互に壁面１０Ｆによって隔てられていないと判定されてしまう。以下に説明する第２実施例では、計算粒子２０ｉと計算粒子２０ｊとが、相互に壁面１０Ｆを隔てて配置されていると判定されるように判定処理がより好適化される。

【0060】

図７は、第２実施例によるシミュレーション方法で用いられる極小距離ベクトルｄ_２ｉについて説明するための模式図である。解析空間に、直交等間隔格子の複数の格子点を定義する。着目する格子点１５ｉの位置ベクトルをｘ_ｉと標記する。着目する格子点の参照範囲２３ｉ内の複数の格子点１５ｊのそれぞれの位置ベクトルをｘ_ｊと標記する。参照範囲２３ｉは、格子点１５ｉの位置を中心とする半径Ｒの円である。半径Ｒは、例えば格子間隔の１０倍程度の大きさとすることができる。

【0061】

格子点１５ｉ、１５ｊのそれぞれの最短距離ベクトルｄ_ｉ、ｄ_ｊが、符号付距離関数から求められる。図７では、格子点１５ｉの最短距離ベクトルｄ_ｉの始点は、壁面１０Ｆの第１領域Ａ１上に位置し、格子点１５ｊの最短距離ベクトルｄ_ｊの始点は、壁面１０Ｆの第２領域Ａ２上に位置する例が示されている。例えば、第１領域Ａ１と第２領域Ａ２とは、相互に直交している。

【0062】

まず、格子点１５ｊに対して次のベクトルｄ_ｉｊを定義する。

【数27】

【0063】

以下に示す評価値Ｄ_ｉｊを計算する。

【数28】

ここで、ベクトルｎ_ｉ、ｎ_ｊ、ｎ_ｉｊは、それぞれベクトルｄ_ｉ、ｄ_ｊ、ｄ_ｉｊと同一の向きを持つ単位ベクトルである。参照範囲２３ｉ内のすべての格子点１５ｊについて評価値Ｄ_ｉｊを計算する。

【0064】

計算結果に基づいて、最も大きな評価値Ｄ_ｉｊを与える格子点１５ｊにおけるベクトルｄ_ｉｊを、着目する格子点１５ｉに対する壁面１５Ｆからの極小距離ベクトルｄ_２ｉとみなす。最短距離ベクトルｄ_ｉは、壁面１５Ｆのうち、格子点１５ｉに最も近い平面状の第１領域Ａ１までの最短距離及び最短経路の向きに相当する。極小距離ベクトルｄ_２ｉは、壁面１５Ｆのうち格子点１５ｉから２番目に近い第２領域Ａ２までの経路の長さの極小値に等しい長さ、及びその経路の向きを表していると考えることができる。

【0065】

次に、式（２８）の意味について説明する。
式（２８）の１番目の中かっこ内の式の値は、ベクトルｄ_ｉとベクトルｄ_ｊとが相互に直交するときに最大値をとり、直交関係からずれるにしたがってその値が小さくなる。式（２８）の２番目の中かっこ内の式の値は、ベクトルｄ_ｉｊとベクトルｄ_ｉとが相互に直交するときに最大値をとり、直交関係からずれるにしたがってその値が小さくなる。式（２８）の小かっこ内の式の値は、ベクトルｄ_ｉｊとベクトルｄ_ｊとが相互に平行であるとき、最大値を取る。

【0066】

したがって、評価値Ｄ_ｉｊが最大となる格子点１５ｊを探索することは、ベクトルｄ_ｉとベクトルｄ_ｊとのなす角度が９０°に近く、ベクトルｄ_ｉｊとベクトルｄ_ｉとのなす角度も９０°に近く、ベクトルｄ_ｉｊとベクトルｄ_ｊとのなす角度が０°に近く、ベクトルｄ_ｉｊの長さが短くなるような格子点１５ｊを探索していると考えることができる。

【0067】

ベクトルｄ_ｉとｄ_ｊとのなす角度が９０°に近いという条件から、ベクトルｄ_ｉとベクトルｄ_ｊとが平行になるような格子点１５ｊは探索対象から除外される。すなわち、壁面１５Ｆの第１領域Ａ１が、最も近い壁面であるような格子点１５ｊが、探索対象から除外される。ベクトルｄ_ｉｊとｄ_ｉとのなす角度も９０°に近いという条件、及びベクトルｄ_ｉｊとｄ_ｊとのなす角度が０°に近いという条件は、着目する格子点１５ｉを通り、最短距離ベクトルｄ_ｉに対して直交する平面に近い位置の格子点１５ｊを探索することと考えられる。したがって、これらの条件を満たす格子点１５ｊに対するベクトルｄ_ｉｊを、着目する格子点１５ｉの極小距離ベクトルｄ_２ｉと見なすことができる。

【0068】

解析空間内の全ての格子点のそれぞれを着目する格子点１５ｉとして、上記操作を行うことにより、すべての格子点について、最短距離ベクトルｄ_ｉ及び極小距離ベクトルｄ_２ｉを求めることができる。

【0069】

次に、図８を参照して、２つの計算粒子が、相互に薄い壁を隔てて配置されているか否かを判定する方法について説明する。図８は、２つの計算粒子２０ｉ、２０ｊが、外側の壁１５Ａの壁面１５Ｆの近くに、かつ薄い壁１０Ｂを隔てて配置されている場合の模式図である。外側の壁１５Ａの壁面１５Ｆを、壁面の第１領域Ａ１といい、薄い壁１５Ｂの壁面１５Ｆを、壁面の第２領域Ａ２ということとする。

【0070】

まず、計算粒子２０ｉ、２０ｊの位置における最短距離ベクトルｄ_ｉ、ｄ_ｊ、及び極小距離ベクトルｄ_２ｉ、ｄ_２ｊを求める。これらのベクトルは、計算粒子２０ｉ、２０ｊの近傍の複数の格子点について計算された最短距離ベクトル及び極小距離ベクトルに基づいて、線形補間演算を行うことによる計算することができる。

【0071】

次に、計算粒子２０ｉに対する壁面の位置ベクトルｘ_ｗｉ、ｘ_ｗ２ｉを、以下の式により計算する。位置ベクトルｘ_ｗｉ及びｘ_ｗ２ｉで表される点は、それぞれ第１領域Ａ１及び第２領域Ａ２上に位置する。

【数29】

【数30】

計算粒子２０ｊについても同様に、位置ベクトルｘ_ｗｊ、ｘ_ｗ２ｊを計算する。位置ベクトルｘ_ｗｊ及びｘ_ｗ２ｊで表される点は、それぞれ第１領域Ａ１及び第２領域Ａ２上に位置する。

【0072】

式（２５）、式（２６）に加えて、以下の２つの式（３１）、式（３２）の合計４つの式のいずれかが満たされる場合、計算粒子２０ｉと計算粒子２０ｊとは、壁面を隔てて配置されていると判定する。

【数31】

【数32】

【0073】

式（３１）は、ベクトルｘ_ｉ－ｘ_ｗ２ｊとベクトルｄ_２ｊとのなす角度θ_２ｊｉが、所定の角度以上であり、式（３２）は、ベクトルｘ_ｊ－ｘ_ｗ２ｉとベクトルｄ_２ｉとのなす角度θ_２ｉｊが、所定の角度以上であることを意味する。パラメータγの値は、式（２５）、式（２６）のパラメータγと同一の値にするとよい。

【0074】

図８に示した例で、式（２５）、式（２６）のみを用いた判定では、薄い壁１５Ｂが検出されず、２つの計算粒子２０ｉ、２０ｊが壁面によって隔てられていないと判定されてしまう。これに対して第２実施例では、式（３１）及び式（３２）が満たされるため、２つの計算粒子２０ｉ、２０ｊは、相互に壁１５Ｂを隔てて配置されていると判定される。

【0075】

図９は、第２実施例によるシミュレーション方法の手順を示すフローチャートである。以下、図４に示したフローチャートと異なるステップについて説明する。図４に示したステップＳ２では、符号付距離関数を算出するが、図８に示したステップＳ２ａでは、符号付距離関数を算出し、さらに、図７を参照して説明した方法で、複数の格子点のそれぞれについて最短距離ベクトルｄ_ｉ及び極小距離ベクトルｄ_２ｉを算出する。

【0076】

図４に示したステップＳ５及びステップＳ９では、着目する計算粒子に対して壁面を隔てて配置された他の計算粒子を近接粒子リストから除外する際に、式（２５）及び式（２６）を用いる。これに対して図８に示したステップＳ５ａ及びステップＳ９ａでは、式（２５）及び式（２６）に加えて式（３１）及び（３２）を用いる。

【0077】

次に、第２実施例の優れた効果について説明する。
第２実施例では、図８に示したように壁面１５Ｆに角部が存在する場合でも、２つの計算粒子が壁面を隔てて配置されているか否かを正しく判定することができる。

【0078】

次に、図１０Ａ～図１０Ｄを参照して、第２実施例で用いられる２つの計算粒子が壁面を隔てて配置されているか否かを判定する方法で実際に判定を行った結果について説明する。図１０Ａ～図１０Ｄは、着目する計算粒子２０ｉに対して参照すべき他の計算粒子２０ｊを示す図である。流体が収容される２つの空間１１Ａ、１１Ｂが、薄い壁１５Ｂで相互に隔てられている。

【0079】

図１０Ａ、図１０Ｂは、一方の空間１１Ａ内の計算粒子２０ｉに着目し、図１０Ｃ、図１０Ｄは、他方の空間１１Ｂ内の計算粒子２０ｉに着目した場合の判定結果を示す。また、図１０Ａ、図１０Ｃは、着目する計算粒子２０ｉが、薄い壁１５Ｂの中央近傍に位置し、図１０Ｂ、図１０Ｄは、着目する計算粒子２０ｉが角部の近傍に位置する場合の判定結果を示す。

【0080】

相対的に薄い灰色の計算粒子２０ｊが、参照すべき他の計算粒子２０ｊとして近接粒子リストから除外されなかったものを示している。なお、参照範囲２３ｉの半径Ｒを、計算粒子の直径の１０倍に設定した。着目する計算粒子２０ｉに対して薄い壁１５Ｂを隔てて配置されている参照範囲２３ｉ内の計算粒子（相対的に濃い色で表されている）が、近接粒子リストから除外されていることが確認された。

【0081】

次に、第２実施例の変形例について説明する。
第２実施例（図７、図８）では、解析空間が二次元平面であり、着目する格子点１５ｉに対して１つの最短距離ベクトルｄ_ｉ及び１つの極小距離ベクトルｄ_２ｉが求まる。解析空間が三次元の場合は、例えば直方体の頂点の近傍のように、壁面の３つの平面の領域で角部が形成される場合がある。この場合は、着目する格子点に対して、最短距離ベクトルｄ_ｉ及び極小距離ベクトルｄ_２ｉの他に、３番目に短い極小距離ベクトルｄ_３ｉが定義される場合がある。

【0082】

解析空間を三次元に拡張する場合は、式（２８）に加えて、以下の式によりＤ２_ｉｊを計算する。

【数33】

ここで、ベクトルｎ_２ｉは、極小距離ベクトルｄ_２ｉと同じ向きの単位ベクトルである。

【0083】

最も大きな評価値Ｄ２_ｉｊを与える格子点１５ｊにおけるベクトルｄ_ｉｊを、着目する格子点１５ｉに対する壁面１５Ｆからの３番目に短い極小距離ベクトルｄ_３ｉとみなすことができる。

【0084】

式（２９）及び式（３０）に加えて、以下の式により位置ベクトルｘ_ｗ３ｉを計算する。

【数34】

【0085】

式（２５）、式（２６）、式（３１）、式（３２）に加えて、以下の２つの式（３５）、式（３６）の合計６個の式のいずれかが満たされる場合、計算粒子２０ｉと計算粒子２０ｊとは、壁面を隔てて配置されていると判定する。

【数35】

【数36】

パラメータγの値は、式（２５）、式（２６）、式（３１）、式（３２）のパラメータγと同一の値にするとよい。

【0086】

上述の判定方法を用いることにより、三次元の流れ場を高精度に解析することが可能になる。

【0087】

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

【符号の説明】

【0088】

１０ 物体
１０Ａ 外側の壁
１０Ｂ 薄い壁
１０Ｆ 壁面
１１Ａ、１１Ｂ 薄い壁で隔てられた空間
１５ｉ 着目する格子点
１５ｊ 着目する格子点の参照範囲内の格子点
２０ 計算粒子
２０ａ 着目する計算粒子に対して壁面を隔てていない位置に配置された計算粒子
２０ｂ 着目する計算粒子に対して壁面を隔てて配置された計算粒子
２０ｉ 着目する計算粒子
２０ｊ 他の計算粒子
２３ｉ 着目する計算粒子に対する参照範囲
３０ 入力部
３１ 処理部
３２ 記憶部
３３ 出力部

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】