JP6897477B2 - 流体シミュレーションプログラム、流体シミュレーション方法および流体シミュレーション装置 - Google Patents

Description

本発明の実施形態は、流体シミュレーションプログラム、流体シミュレーション方法および流体シミュレーション装置に関する。

従来、流体や弾性体などの連続体（以下、流体）の運動を数値的に計算する流体解析手法の一つとして、連続体を粒子として離散化し、粒子分布で表現する粒子法がある。この粒子法としては、Smoothed Particle Hydrodynamics（SPH）法や Moving Particle Semi-implicitもしくはMoving Particle Simulation（MPS）法が知られている。

これら、粒子分布を用いて解析対象を表現する手法においては、ある粒子から影響半径と呼ばれる距離（ｈ）以内にある粒子を近傍粒子としてその粒子の情報を用いる。一例として、ＳＰＨ法で非圧縮流体の運動方程式として以下の式（１）を解くことを考える。

ここで、左辺は流体要素の加速度、右辺第一項は圧力勾配項、第二項は外力項である。これをＳＰＨ法を用いて離散化する。ｘ_ｉ，ｖ_ｉ，ρ_ｉ，ｍ_ｉ，ｐ_ｉをそれぞれ、粒子ｉの位置ベクトル、速度ベクトル、密度、質量、圧力とすると、時間発展について以下の式（２ａ）、（２ｂ）のとおりに計算することができる。

ここで、Δｔはシミュレーションの時間刻み、上付き添え字ｎは時間ステップ数である。また、角括弧＜＞はＳＰＨ法による重ね合わせ評価を表し、次の式（３）のとおりである。

ここで、和の対象ｊは粒子ｉから半径ｈ以内に存在する粒子(近傍粒子)のインデックスを渡るものとする。また、以後ＳＰＨ法の重ね合わせの和をとるときはすべて同様である。Ｗはカーネル関数で、粒子の分布から連続場を構成するのに用い、以下の式（４）の３次のスプライン関数等がよく使われる。ここで、βは規格化係数であり、解析対象が三次元の場合はπh³ である。

各粒子の圧力ｐ_ｉは非圧縮条件∇・ｖ_ｉ＝０を満たすように決めるため、式（５）の圧力ポアソン方程式を陰的に解くことで求める。

ここで、ρ_０は物質の基準密度、ｖ_ｉ ^＊は外力項によって更新された中間状態の速度であり、次の式（６）で表現される。

式（５）の左辺は圧力のラプラシアンであり、次の式（７）で離散化される。

また、式（５）の右辺の速度発散は、次の式（８）で離散化される。

このように離散化することで、式（７）の圧力ポアソン方程式は、次の式（９）、（１０）として表すことができ、次の式（１１）という線形方程式に帰着する。

上記の圧力ポアソン方程式は、ある粒子の近傍に十分な粒子がない時に、その粒子を自由表面とみなし、自由表面の粒子に対して圧力を０と定めて適切な境界条件を与えることで、共役勾配法等の反復法により解くことができる。

特開２０１０−２４３２９３号公報 特開２０１６−１２５９３４号公報

J.J. Monaghan, "Smoothed Particle Hydrodynamics", Annu. Rev. Astron. Astrophys. 30: 543-74(1992) 浅井光輝，藤本啓介，田邊将一，別府万寿博，"階段状の非適合境界を有する粒子法解析における仮想マーカーを用いたすべり・非すべり境界処理法"，Transactions of JSCES, Paper No.20130011，２０１３年

しかしながら、上記の従来技術では、細い隙間に流体が入り込むことで適切な境界条件が設定されない場合に、解が一意に定まらず、圧力が発散するなどの不具合が生じることがある。例えば、細い隙間に流体が入り込んだ時には、全ての流体粒子に対して近傍粒子が存在することになり、狭い隙間において自由表面と判定できる粒子が存在しなくなることから、適切な境界条件が設定されないこととなる。

１つの側面では、粒子法における解の収束性を確保することを可能とする流体シミュレーションプログラム、流体シミュレーション方法および流体シミュレーション装置を提供することを目的とする。

第１の案では、流体シミュレーションプログラムは、抽出する処理と、求める処理と、算出する処理とをコンピュータに実行させる。抽出する処理は、所定の時間における各粒子の粒子データに含まれる位置情報をもとに、壁との境界にかかる壁境界粒子の中から、流体粒子との距離が所定値以下の粒子を第１粒子とし、第１粒子との距離が所定値以下の粒子を第２粒子として抽出する。求める処理は、抽出した第１粒子および第２粒子をもとに、各粒子の圧力を求める圧力ポアソン方程式の境界条件を設定して各粒子にかかる圧力を求める。算出する処理は、求めた圧力に基づいて次の時間における各粒子の粒子データを算出する。

本発明の１実施態様によれば、粒子法における解の収束性を確保することができる。

図１は、実施形態にかかる流体シミュレーション装置の構成例を示す機能ブロック図である。 図２−１は、自由流体粒子の入力データの一例を示す説明図である。 図２−２は、壁境界粒子の入力データの一例を示す説明図である。 図２−３は、壁内部粒子の入力データの一例を示す説明図である。 図３−１は、自由流体粒子の中間データの一例を示す説明図である。 図３−２は、壁境界粒子の中間データの一例を示す説明図である。 図３−３は、壁内部粒子の中間データの一例を示す説明図である。 図４は、実施形態にかかる流体シミュレーション装置の動作例を示すフローチャートである。 図５は、近傍リストの一例を示す説明図である。 図６−１は、粒子の近傍を表す空間領域の分類を説明する説明図である。 図６−２は、表面粒子と内部粒子の抽出例を説明する説明図である。 図７は、粒子の分類処理を例示するフローチャートである。 図８は、粒子の分類を説明する説明図である。 図９は、求解処理を例示するフローチャートである。 図１０−１は、解が発散するケースを説明する説明図である。 図１０−２は、境界条件の改善例を説明する説明図である。 図１１は、実施形態にかかる流体シミュレーション装置のハードウエア構成例を示すブロック図である。

以下、図面を参照して、実施形態にかかる流体シミュレーションプログラム、流体シミュレーション方法および流体シミュレーション装置を説明する。実施形態において同一の機能を有する構成には同一の符号を付し、重複する説明は省略する。なお、以下の実施形態で説明する流体シミュレーションプログラム、流体シミュレーション方法および流体シミュレーション装置は、一例を示すに過ぎず、実施形態を限定するものではない。また、以下の各実施形態は、矛盾しない範囲内で適宜組みあわせてもよい。

図１は、実施形態にかかる流体シミュレーション装置の構成例を示す機能ブロック図である。図１に示すように、流体シミュレーション装置１は、粒子法による流体シミュレーションを実行するワークステーション等の情報処理装置であり、制御部１０と、記憶部２０とを有する。

制御部１０は、ＣＰＵ（Central Processing Unit）などの電子回路に対応する。制御部１０は、各種の処理手順を規定したプログラムや制御データを格納するための内部メモリを有し、これらによって種々の処理を実行する。制御部１０は、シミュレーション受付部１１と、粒子分類部１２と、演算部１３と、出力部１４とを有する。

記憶部２０は、例えば、ＲＡＭ、フラッシュメモリ（Flash Memory）などの半導体メモリ素子、または、ハードディスク、光ディスクなどの記憶装置である。記憶部２０は、粒子情報２１と、逐次結果情報２２とを有する。

粒子情報２１は、流体シミュレーションの設定条件として入力される入力データであり、主には流体シミュレーションにかかる各粒子の情報である。この粒子情報２１における入力データは、例えば計算パラメータと、３種類の各粒子の粒子データとを有する。

計算パラメータは、シミュレーション対象となる流体の物性値やシミュレーション環境などを示すパラメータである。例えば、計算パラメータとしては、各粒子が影響を及ぼす影響半径、各粒子に及ぶ外力ベクトル、シミュレーションの時間刻み幅、シミュレーション対象の物質の基準密度、一粒子の質量、ゼロ除算防止パラメータηなどを含む。

粒子データは、自由流体粒子、壁境界粒子、壁内部粒子の３種類を含む。図２−１は、自由流体粒子の入力データの一例を示す説明図である。図２−２は、壁境界粒子の入力データの一例を示す説明図である。図２−３は、壁内部粒子の入力データの一例を示す説明図である。

自由流体粒子は、シミュレーションの解析対象となる流体を表す粒子である。すなわち、自由流体粒子は、流体粒子の一例である。図２−１に示すように、自由流体粒子にかかる自由流体粒子入力データ１００は、各粒子（ｉ、ｉ＋１、…）の位置座標１０１、速度１０２密度１０３および質量１０４を有する。なお、位置座標１０１および速度１０２については、入力データを初期値として、時間ステップごとに計算された位置および速度の値で更新される。

壁境界粒子は、自由流体粒子の移動が妨げられる固定境界表面に配置し、圧力ポアソン方程式の境界条件を表すための粒子である。図２−２に示すように、壁境界粒子にかかる壁境界粒子入力データ２００は、各粒子（ｉ、ｉ＋１、…）の位置座標２０１、密度２０２および質量２０３を有する。

壁内部粒子は、固定境界内部に配置し、速度発散と圧力勾配を計算する際に用いる粒子である。図２−３に示すように、壁内部粒子にかかる壁内部粒子入力データ３００は、各粒子（ｉ、ｉ＋１、…）の位置座標３０１、法線ベクトル３０２、固定境界面までの距離３０３、密度３０４および質量３０５を有する。なお、これらの粒子の密度、質量については、計算パラメータに含まれる値を使用してもよい。

逐次結果情報２２は、流体シミュレーションにおいて、各粒子の時間ステップごとのシミュレーション結果を示す中間データである。図３−１は、自由流体粒子の中間データの一例を示す説明図である。図３−２は、壁境界粒子の中間データの一例を示す説明図である。図３−３は、壁内部粒子の中間データの一例を示す説明図である。

図３−１に示すように、自由流体粒子にかかる自由流体粒子中間データ１１０は、各粒子（ｉ、ｉ＋１、…）の圧力１１１、中間速度１１２、速度発散１１３および境界フラグ１１４を有する。ここで、境界フラグ１１４は、各自由流体粒子について、内部粒子または表面粒子の分類を示すフラグである。

図３−２に示すように、壁境界粒子にかかる壁境界粒子中間データ２１０は、各粒子（ｉ、ｉ＋１、…）の速度発散２１１、圧力２１２および境界フラグ２１３を有する。ここで、境界フラグ２１３は、各壁境界粒子について、固定流体粒子、圧力境界粒子または非使用粒子の分類を示すフラグである。

図３−３に示すように、壁内部粒子にかかる壁内部粒子中間データ３１０は、各粒子（ｉ、ｉ＋１、…）の仮想マーカ位置３１１、仮想マーカ速度３１２、仮想マーカ加速度３１３、仮想マーカ圧力３１４、中間速度３１５および圧力３１６を有する。ここで、仮想マーカ位置３１１は、壁内部粒子の固定境界面に対する鏡像位置であり、仮想マーカ速度３１２、仮想マーカ加速度３１３および仮想マーカ圧力３１４は、鏡像位置における速度、加速度および圧力の仮想値である。

シミュレーション受付部１１は、流体シミュレーションの実行要求を受け付ける。例えば、シミュレーション受付部１１は、計算パラメータ、自由流体粒子入力データ１００、壁境界粒子入力データ２００および壁内部粒子入力データ３００などの粒子情報２１を含む実行要求を受け付ける。次いで、シミュレーション受付部１１は、受け付けた粒子情報２１を記憶部２０に格納し、シミュレーション受付部１１、粒子分類部１２および出力部１４による流体シミュレーションの処理を開始する。

粒子分類部１２は、計算パラメータの時間刻み幅による時間ステップごとに、粒子情報２１における各粒子の位置情報（位置座標１０１、２０１、３０１）をもとに、各粒子の分類を行う。具体的には、粒子分類部１２は、自由流体粒子について、表面粒子と、内部粒子とに分類する。また、粒子分類部１２は、壁境界粒子について、固定流体粒子と、圧力境界粒子と、非使用粒子とに分類する。

演算部１３は、粒子分類部１２による各粒子の分類結果をもとに、各粒子の圧力を求める圧力ポアソン方程式の境界条件を設定し、各粒子にかかる圧力を求める。具体的には、演算部１３は、壁境界粒子について、固定流体粒子として分類（抽出）された粒子については速度を所定値として固定する流体粒子とし、圧力ポアソン方程式を解く際の自由度に含めるようにする。また、演算部１３は、壁境界粒子について、圧力境界粒子として分類（抽出）された粒子については圧力を０とする境界条件を設定する。次いで、演算部１３は、設定した条件で圧力ポアソン方程式により各粒子にかかる圧力を算出する。次いで、演算部１３は、求めた圧力に基づいて次の時間における各粒子の粒子データ（例えば速度および位置）を算出する。

出力部１４は、演算部１３の算出結果を出力する。具体的には、出力部１４は、時間ステップごとに算出された算出結果を逐次結果情報２２として記憶部２０に格納する。なお、出力部１４は、シミュレーションの開始から終了までの時間ステップごとに算出された逐次結果情報２２に基づくシミュレーション結果を、ディスプレイなどの表示装置への表示出力やファイル出力してもよい。

ここで、流体シミュレーション装置１における流体シミュレーションの処理内容について、詳細を説明する。図４は、実施形態にかかる流体シミュレーション装置１の動作例を示すフローチャートである。

図４に示すように、処理が開始されると、流体シミュレーション装置１のシミュレーション受付部１１は、ユーザからの入力操作などを介して、粒子情報２１にかかる入力データの取得を行う（Ｓ１）。

次いで、流体シミュレーション装置１の粒子分類部１２、演算部１３および出力部１４は、粒子情報２１および逐次結果情報２２を参照し、シミュレーションの終了時間に対応するステップ（時間ステップ）数分のループ処理（Ｓ２〜Ｓ１２）を実行する。

Ｓ２に次いで、粒子分類部１２は、粒子情報２１における位置情報（位置座標１０１、２０１、３０１）に基づいて、粒子ごとの、近傍にある粒子の近傍リストを作成する（Ｓ３）。具体的には、粒子分類部１２は、各粒子において、粒子の位置から計算パラメータにおける影響半径内の距離にある粒子を抽出し、近傍の粒子を識別する粒子番号の近傍リストを作成する。

図５は、近傍リストの一例を示す説明図である。図５に示すように、近傍リスト３０は、粒子（…ｉ−１、ｉ、ｉ＋１、ｉ＋２、…）ごとに、影響半径内にある近傍の粒子を並べたものである。例えば、粒子（ｉ）については、粒子（ｊ_１、ｊ_２、…）が近傍にある粒子のリストとして挙げられている。また、粒子（ｉ＋１）については、粒子（ｊ’_１，ｊ’_２、…）が近傍にある粒子のリストとして挙げられている。

次いで、演算部１３は、自由流体粒子の各粒子の速度（ｖ_ｉ ^ｎ）に対して、外力ベクトル（ｇ）と時間刻み幅（Δｔ）とによる速度更新（式（６））を行い、中間速度（ｖ_ｉ ^＊）を算出する（Ｓ４）。

次いで、演算部１３は、自由流体粒子に対して、速度発散を算出する（Ｓ５）。これは、仮想マーカ位置３１１における速度評価、壁内部粒子の速度評価、速度発散の算出、という手順で行う。

具体的には、演算部１３は、仮想マーカ位置３１１を壁内部粒子の固定境界面に対する鏡像位置として算出する。例えば、ある壁内部粒子の仮想マーカ位置３１１は、粒子の位置座標３０１、法線ベクトル３０２および固定境界面までの距離３０３を用いて次の式（１２）で算出される。

この位置における速度（仮想マーカ速度３１２）は、次の式（１３）によって評価する。ここで、和をとる範囲は自由流体粒子の近傍粒子とする。

演算部１３は、壁内部粒子の中間速度３１５を、仮想マーカ位置３１１の仮想マーカ速度３１２からｖ_ｉ ^＊＝−ｖ_ｉ ^ＶＭとして算出する。これは、固定境界面において非すべり条件を課したことに相当する。

次いで、演算部１３は、各自由流体粒子に対して速度発散を式（８）用いて算出する。ここで、和をとる範囲は、自由流体粒子と壁内部粒子の近傍粒子の範囲とする。

次いで、粒子分類部１２は、自由流体粒子の中から表面粒子とする粒子を検出（抽出）する（Ｓ６）。具体的には、粒子分類部１２は、ある粒子に注目したとき、いずれかの方向に近傍粒子（近傍リスト３０にある、その粒子からの距離が影響半径よりも小さい粒子）が存在しない時を表面粒子とする。また、粒子分類部１２は、全ての方向に近傍粒子が存在するときを内部粒子とする。

図６−１は、粒子の近傍を表す空間領域の分類を説明する説明図である。図６−２は、表面粒子と内部粒子の抽出例を説明する説明図である。

図６−１に示すように、２次元を例にして説明すると、先ず、注目する粒子（ｉ）を中心として影響半径内の空間領域をｘ＋、ｘ−、ｙ＋、ｙ−の４領域に分類する。次いで、粒子分類部１２は、影響半径内の粒子がいずれの空間領域に属するかを判断し、粒子が一つも所属しない領域があった場合、粒子（ｉ）を表面粒子であると抽出する。

図６−２に示すように、ケースＣ１では、ｙ＋の領域に近傍粒子が存在しないため、中心の粒子（ｉａ）については表面粒子として抽出される。一方、ケースＣ２では、全ての領域に近傍粒子が存在するため、粒子（ｉｂ）は内部粒子として抽出される。粒子分類部１２は、この判定結果について、自由流体粒子中間データ１１０の境界フラグ１１４に分類結果として保存する。

次いで、粒子分類部１２は、粒子情報２１における位置座標１０１、２０１、３０１に基づく粒子の分類処理により、壁境界粒子を固定流体粒子、圧力境界粒子、非使用粒子に分類する（Ｓ７）。

図７は、粒子の分類処理を例示するフローチャートである。図７に示すように、分類処理が開始されると、粒子分類部１２は、全ての壁境界粒子の境界フラグ２１３を非使用粒子とする（Ｓ２０）。

次いで、粒子分類部１２は、自由流体粒子ごとのループ処理（ｉ）を実行する（Ｓ２１〜Ｓ２８）。ループ処理（ｉ）が開始されると、粒子分類部１２は、自由流体粒子を一つ選び、粒子ｉとする（Ｓ２２）。

次いで、粒子分類部１２は、壁境界粒子ごとのループ処理（ｊ）を実行する（Ｓ２３〜Ｓ２７）。ループ処理（ｊ）が開始されると、粒子分類部１２は、壁境界粒子を一つ選び粒子ｊとする（Ｓ２４）。

次いで、粒子分類部１２は、粒子ｉの位置座標１０１と、粒子ｊの位置座標２０１とを用いて粒子ｉと粒子ｊの距離（ｒ）を求め、粒子ｉと粒子ｊの距離が影響半径以下であるか否かを判定する（Ｓ２５）。

影響半径以下である場合（Ｓ２５：ＹＥＳ）、粒子分類部１２は、粒子ｊの境界フラグ２１３を固定流体粒子とする（Ｓ２６）。影響半径以下でない場合（Ｓ２５：ＮＯ）、粒子分類部１２は、Ｓ２６をスキップして、境界フラグ２１３をそのままとする。

次いで、粒子分類部１２は、まだ選択されていない壁境界粒子があれば、その粒子を粒子ｊとして、ループ処理（ｊ）を繰り返す。また、粒子分類部１２は、ループ処理（ｊ）の後、まだ選択されていない自由流体粒子があれば、その粒子を粒子ｉとして、ループ処理（ｉ）を繰り返す。

次いで、粒子分類部１２は、固定流体粒子ごとのループ処理（ｊ’）を実行する（Ｓ２９〜Ｓ３７）。ループ処理（ｊ’）が開始されると、粒子分類部１２は、境界フラグ２１３を固定流体粒子とされた壁境界粒子を一つ選び粒子ｊ’とする（Ｓ３０）。

次いで、粒子分類部１２は、壁境界粒子ごとのループ処理（ｊ’’）を実行する（Ｓ３１〜Ｓ３６）。ループ処理（ｊ’’）が開始されると、粒子分類部１２は、壁境界粒子を一つ選び粒子ｊ’’とする（Ｓ３２）。

次いで、粒子分類部１２は、粒子ｊ’’の境界フラグ２１３が非使用粒子であるか否かを判定する（Ｓ３３）。非使用粒子でない場合（Ｓ３３：ＮＯ）、粒子分類部１２は、Ｓ３６へ処理をスキップし、次の粒子ｊ’’を選ぶループ処理を行う。

非使用粒子である場合（Ｓ３３：ＹＥＳ）、粒子分類部１２は、粒子ｊ’、ｊ’’の位置座標２０１を用いて粒子ｊ’、ｊ’’の距離（ｒ_{ｊ’ｊ’’}）を求め、粒子ｉと粒子ｊの距離が影響半径以下であるか否かを判定する（Ｓ３４）。

影響半径以下である場合（Ｓ３４：ＹＥＳ）、粒子分類部１２は、粒子ｊ’’の境界フラグ２１３を圧力境界粒子とする（Ｓ３５）。影響半径以下でない場合（Ｓ３４：ＮＯ）、粒子分類部１２は、Ｓ３５をスキップして、境界フラグ２１３をそのままとする。

次いで、粒子分類部１２は、まだ選択されていない壁境界粒子があれば、その粒子を粒子ｊ’’として、ループ処理（ｊ’’）を繰り返す。また、粒子分類部１２は、ループ処理（ｊ’’）の後、まだ選択されていない固定流体粒子があれば、その粒子を粒子ｊ’として、ループ処理（ｊ’）を繰り返す。

図８は、粒子の分類を説明する説明図である。図８に示すように自由流体粒子４０については、Ｓ６の処理により、内部粒子４１または表面粒子４２に分類される。また、壁境界粒子５０については、Ｓ７の分類処理により、固定流体粒子５１、圧力境界粒子５２および非使用粒子５３のいずれかに分類される。

図４に戻り、Ｓ７に次ぐ処理を説明する。なお、非使用粒子５３については、このステップが終わるまで使用されないので、以後、無視することとする。ここで改めて自由流体粒子４０、固定流体粒子５１、圧力境界粒子５２などに通し番号を振り、以後、その番号を用いて粒子ｉ、粒子ｊ等として参照する。また、演算部１３は、固定流体粒子５１、圧力境界粒子５２の圧力２１２および速度発散２１１を０とする。

Ｓ７に次いで、演算部１３は、各粒子の圧力を求める圧力ポアソン方程式として、式（１１）に示した線形方程式を解くことで各粒子の圧力を求める（Ｓ８）。

具体的には、演算部１３は、作業データとして係数行列Ａと右辺ベクトルｂを用意する。係数行列Ａは粒子数の２乗個の実数値が格納できれば十分であるが、そのほとんどが０値となるので、非ゼロ要素のみを効率よく格納するためにCompressed Sparse Row（CSR）形式などを用いてもよい。右辺ベクトルｂは粒子数と同数の実数値を格納できる配列として確保する。演算部１３は、この係数行列Ａと右辺ベクトルｂに対して行列方程式（式１１）を解く求解処理を行う。

図９は、求解処理を例示するフローチャートである。図９に示すように、処理が開始されると、演算部１３は、係数行列Ａ_ｉｊの値を全て０とし、全要素０で初期化する（Ｓ４０）。

次いで、演算部１３は、自由流体粒子または壁境界粒子を一つ選び粒子ｉとする（Ｓ４１）。次いで、演算部１３は、粒子ｉの境界フラグ１１４または境界フラグ２１３が内部粒子４１または固定流体粒子５１であるか否かを判定する（Ｓ４２）。

内部粒子４１または固定流体粒子５１でない場合（Ｓ４２：ＮＯ）、演算部１３は、係数行列Ａ_ｉｊの要素を１とする（Ｓ４９）。次いで、演算部１３は、右辺ベクトルｂ_ｉの要素の値を粒子ｉの圧力ｐ_ｉ（圧力１１１または圧力２１２）とし（Ｓ５０）、Ｓ５１へ処理を進める。

内部粒子４１または固定流体粒子５１である場合（Ｓ４２：ＹＥＳ）、演算部１３は、粒子ｉの近傍リスト３０に含まれる粒子ごとのループ処理（ｊ）を実行する（Ｓ４３〜Ｓ４７）。ループ処理（ｊ）が開始されると、演算部１３は、粒子ｉの近傍リスト３０に含まれる粒子を一つ選び、粒子ｊとする（Ｓ４４）。

次いで、演算部１３は、係数行列Ａ_ｉｊの要素について、式（９）の右辺上によって値を算出する（Ｓ４５）。ここで、演算部１３は、各粒子の位置座標１０１または位置座標２０１、密度１０３または密度２０２、質量１０４または質量２０３を用いる。次いで、演算部１３は、係数行列Ａ_ｉｉの要素に、Ｓ４５で算出した値を加算する（Ｓ４６）。

次いで、演算部１３は、粒子ｉの近傍リスト３０に含まれる粒子の中でまだ選択されていない粒子があれば、その粒子を粒子ｊとして、ループ処理（ｊ）を繰り返す。

粒子ｉの近傍粒子の全てに対してＳ４３〜Ｓ４７の処理が済んだ後、演算部１３は、右辺ベクトルｂ_ｉの値を式（１０）を用いて算出する（Ｓ４８）。ここでは、各粒子の速度発散１１３または速度発散２１１を用いる。

次いで、演算部１３は、まだ選んでいない内部粒子４１か固定流体粒子５１があるか否かを判定する（Ｓ５１）。まだ選んでいない内部粒子４１か固定流体粒子５１がある場合（Ｓ５１：ＹＥＳ）、演算部１３は、Ｓ４１へ処理を戻し、新たな粒子ｉを選ぶ。

まだ選んでいない内部粒子４１か固定流体粒子５１がない場合（Ｓ５１：ＮＯ）、演算部１３は、ここまでの処理で得られた係数行列Ａ_ｉｊと右辺ベクトルｂおよび各粒子の圧力１１１または圧力２１２をｐとし、行列方程式を解く（Ｓ５２）。具体的には、演算部１３は、行列方程式Ａ_ｉｊｐ_ｊ＝ｂ_ｉを解くことで、圧力を求める。この求解には、例えば共役勾配法などの解法を用いればよい。

図４に戻り、Ｓ８に次いで、演算部１３は、Ｓ８の求解処理による算出結果をもとに、各粒子について圧力勾配による加速度（速度）の算出を行う（Ｓ９）。具体的には、演算部１３は、Ｓ８により各自由流体粒子の圧力ｐ_ｉが求められたところで、次に壁内部粒子の圧力を求める。

演算部１３は、壁内部粒子の圧力を求めるために、先ず仮想マーカ加速度３１３と仮想マーカ圧力３１４を次の式（１４）、（１５）によって求める。

次に、演算部１３は、式（１４）、（１５）によって求めた値を用いて、壁内部粒子の圧力ｐ_ｉ ^ｗａｌｌを次の式（１６）によって算出する。

このようにして自由流体粒子と壁内部粒子の圧力が求められたところで、演算部１３は、式（３）を用いて圧力勾配を求め、圧力勾配により粒子の加速度、速度を更新する。なお、式（３）の和については、自由流体粒子と壁内部粒子との近傍をわたるものとする。

次いで、演算部１３は、速度・加速度の時間積分を行い、各粒子の位置を更新して時間を１ステップ進める（Ｓ１０）。具体的には、演算部１３は、求めた速度に１ステップ分の時間をかけ合わせて各粒子の移動距離を求め、各粒子の位置（位置座標１０１）を更新することで１ステップ分の計算が完了する。

次いで、出力部１４は、計算結果を逐次結果情報２２としてファイル出力する（Ｓ１１）。流体シミュレーション装置１は、このＳ２〜Ｓ１２のループ処理を繰り返すことで、シミュレーションの開始から終了までの、計算パラメータの時間刻み幅による時間ステップごとのシミュレーション結果を得る。

以上のように、流体シミュレーション装置１は、粒子分類部１２と、演算部１３とを有する。粒子分類部１２は、各粒子の時間ステップごとの粒子データに含まれる位置情報をもとに、壁境界粒子５０の中から、自由流体粒子４０との距離が所定値以下の粒子を固定流体粒子５１として抽出する。また、粒子分類部１２は、固定流体粒子５１との距離が所定値以下の粒子を圧力境界粒子５２として抽出する。演算部１３は、抽出した固定流体粒子５１および圧力境界粒子５２をもとに、各粒子の圧力を求める圧力ポアソン方程式の境界条件を設定して各粒子にかかる圧力を求める。次いで、演算部１３は、求めた圧力に基づいて次の時間ステップにおける各粒子の粒子データを算出する。このように壁境界粒子５０より抽出した固定流体粒子５１および圧力境界粒子５２をもとに各粒子の圧力を求める圧力ポアソン方程式の境界条件を設定することで、流体シミュレーション装置１は、例えば細い隙間に流体が入り込む場合であっても適切な境界条件を設定することができる。このため、流体シミュレーション装置１は、粒子法における解の収束性を確保することができる。

ここで、適切な境界条件を設定し、境界条件を改善することで、粒子法における解の収束性を確保されることを、解が発散するケースと比較して説明する。図１０−１は、解が発散するケースを説明する説明図である。図１０−２は、境界条件の改善例を説明する説明図である。

図１０−１に示すように、壁境界粒子５０が一律に設定されるケースでは、狭い隙間に入り込む自由流体粒子４０について、粒子４０ａを取り囲む粒子が多く、内部粒子として判定される。このケースでは、粒子４０ａの近傍にある壁境界粒子５０において圧力ポアソン方程式を解く際の境界条件に含める粒子が存在しないことから、境界条件が設定されずに圧力ポアソン方程式を解くことで解が発散する場合がある。

これに対し、本実施形態では、壁境界粒子５０の中から、自由流体粒子４０との距離が所定値以下の粒子を固定流体粒子５１として抽出し、固定流体粒子５１との距離が所定値以下の粒子を圧力境界粒子５２として抽出する。そして、壁境界粒子５０より抽出した固定流体粒子５１および圧力境界粒子５２をもとに各粒子の圧力を求める圧力ポアソン方程式の境界条件を設定する。

例えば、流体シミュレーション装置１は、固定流体粒子５１を速度固定の流体粒子とし、圧力境界粒子５２を圧力を０として圧力ポアソン方程式の境界条件を設定する。このように、流体シミュレーション装置１は、流体の近傍となった固定流体粒子５１を圧力ポアソン方程式を解く際の自由度に含め、固定流体粒子５１近傍の圧力境界粒子５２は圧力０の表面粒子として圧力ポアソン方程式の境界条件とする。これにより、流体シミュレーション装置１は、例えば細い隙間に流体が入る事象を破綻なく計算できる。

具体的には、図１０−２に示すように、自由流体粒子４０の近傍にある壁境界粒子５０については、例えば速度固定の固定流体粒子５１として扱い、圧力ポアソン方程式の自由度に含める。また、流体の近傍ではないが、固定流体粒子５１の近傍にある壁境界粒子を圧力境界粒子５２（表面粒子）として、圧力０と定めることで、圧力ポアソン方程式の境界条件とする。なお、流体の近傍でも固定流体粒子５１の近傍でもない壁境界粒子５０については、非使用粒子５３として圧力ポアソン方程式の算出時に使用しないものとする。このように、固定流体粒子５１および圧力境界粒子５２をもとに各粒子の圧力を求める圧力ポアソン方程式の境界条件を設定することで、狭い隙間に入り込む自由流体粒子４０がある場合でも、解の収束性を確保することができる。

なお、図示した各装置の各構成要素は、必ずしも物理的に図示の如く構成されていることを要しない。すなわち、各装置の分散・統合の具体的形態は図示のものに限られず、その全部または一部を、各種の負荷や使用状況などに応じて、任意の単位で機能的または物理的に分散・統合して構成することができる。

また、流体シミュレーション装置１で行われる各種処理機能は、ＣＰＵ（またはＭＰＵ、ＭＣＵ（Micro Controller Unit）等のマイクロ・コンピュータ）上で、その全部または任意の一部を実行するようにしてもよい。また、各種処理機能は、ＣＰＵ（またはＭＰＵ、ＭＣＵ等のマイクロ・コンピュータ）で解析実行されるプログラム上、またはワイヤードロジックによるハードウエア上で、その全部または任意の一部を実行するようにしてもよいことは言うまでもない。また、流体シミュレーション装置１で行われる各種処理機能は、クラウドコンピューティングにより、複数のコンピュータが協働して実行してもよい。

ところで、上記の実施形態で説明した各種の処理は、予め用意されたプログラムをコンピュータで実行することで実現できる。そこで、以下では、上記の実施形態と同様の機能を有するプログラムを実行するコンピュータ（ハードウエア）の一例を説明する。図１１は、実施形態にかかる流体シミュレーション装置１のハードウエア構成例を示すブロック図である。

図１１に示すように、流体シミュレーション装置１は、各種演算処理を実行するＣＰＵ４０１と、データ入力を受け付ける入力装置４０２と、モニタ４０３と、音声等の出力装置４０４とを有する。また、流体シミュレーション装置１は、記憶媒体からプログラム等を読み取る媒体読取装置４０５と、各種装置と接続するためのインタフェース装置４０６と、有線または無線により外部機器と通信接続するための通信装置４０７とを有する。また、流体シミュレーション装置１は、各種情報を一時記憶するＲＡＭ４０８と、ハードディスク装置４０９とを有する。また、流体シミュレーション装置１内の各部（４０１〜４０９）は、バス４１０に接続される。

ハードディスク装置４０９には、上記の実施形態で説明したシミュレーション受付部１１、粒子分類部１２、演算部１３および出力部１４などで各種の処理を実行するためのプログラム４１１が記憶される。また、ハードディスク装置４０９には、プログラム４１１が参照する各種データ４１２（例えば粒子情報２１、逐次結果情報２２等）が記憶される。入力装置４０２は、例えば、流体シミュレーション装置１の操作者から操作情報の入力を受け付ける。モニタ４０３は、例えば、操作者が操作する各種画面を表示する。インタフェース装置４０６は、例えば印刷装置等が接続される。通信装置４０７は、ＬＡＮ（Local Area Network）等の通信ネットワークと接続され、通信ネットワークを介した外部機器との間で各種情報をやりとりする。

ＣＰＵ４０１は、ハードディスク装置４０９に記憶されたプログラム４１１を読み出して、ＲＡＭ４０８に展開して実行することで、シミュレーション受付部１１、粒子分類部１２、演算部１３および出力部１４にかかる各種の処理を行う。なお、プログラム４１１は、ハードディスク装置４０９に記憶されていなくてもよい。例えば、流体シミュレーション装置１が読み取り可能な記憶媒体に記憶されたプログラム４１１を読み出して実行するようにしてもよい。流体シミュレーション装置１が読み取り可能な記憶媒体は、例えば、ＣＤ−ＲＯＭやＤＶＤディスク、ＵＳＢ（Universal Serial Bus）メモリ等の可搬型記録媒体、フラッシュメモリ等の半導体メモリ、ハードディスクドライブ等が対応する。また、公衆回線、インターネット、ＬＡＮ等に接続された装置にこのプログラム４１１を記憶させておき、流体シミュレーション装置１がこれらからプログラム４１１を読み出して実行するようにしてもよい。

以上の実施形態に関し、さらに以下の付記を開示する。

（付記１）所定の時間における各粒子の粒子データに含まれる位置情報をもとに、壁との境界にかかる壁境界粒子の中から、流体粒子との距離が所定値以下の粒子を第１粒子とし、前記第１粒子との距離が所定値以下の粒子を第２粒子として抽出し、
抽出した前記第１粒子および前記第２粒子をもとに、各粒子の圧力を求める圧力ポアソン方程式の境界条件を設定して各粒子にかかる圧力を求め、
求めた前記圧力に基づいて次の時間における各粒子の粒子データを算出する、
処理をコンピュータに実行させることを特徴とする流体シミュレーションプログラム。

（付記２）前記求める処理は、前記第１粒子を速度を所定値とする流体粒子とし、前記第２粒子の圧力を０として前記境界条件を設定する、
付記１に記載の流体シミュレーションプログラム。

（付記３）所定の時間における各粒子の粒子データに含まれる位置情報をもとに、壁との境界にかかる壁境界粒子の中から、流体粒子との距離が所定値以下の粒子を第１粒子とし、前記第１粒子との距離が所定値以下の粒子を第２粒子として抽出し、
抽出した前記第１粒子および前記第２粒子をもとに、各粒子の圧力を求める圧力ポアソン方程式の境界条件を設定して各粒子にかかる圧力を求め、
求めた前記圧力に基づいて次の時間における各粒子の粒子データを算出する、
処理をコンピュータが実行することを特徴とする流体シミュレーション方法。

（付記４）前記求める処理は、前記第１粒子を速度を所定値とする流体粒子とし、前記第２粒子の圧力を０として前記境界条件を設定する、
付記３に記載の流体シミュレーション方法。

（付記５）所定の時間における各粒子の粒子データに含まれる位置情報をもとに、壁との境界にかかる壁境界粒子の中から、流体粒子との距離が所定値以下の粒子を第１粒子とし、前記第１粒子との距離が所定値以下の粒子を第２粒子として抽出する分類部と、
抽出した前記第１粒子および前記第２粒子をもとに、各粒子の圧力を求める圧力ポアソン方程式の境界条件を設定して各粒子にかかる圧力を求め、求めた前記圧力に基づいて次の時間における各粒子の粒子データを算出する算出部と、
を有することを特徴とする流体シミュレーション装置。

（付記６）前記算出部は、前記第１粒子を速度を所定値とする流体粒子とし、前記第２粒子の圧力を０として前記境界条件を設定する、
付記５に記載の流体シミュレーション装置。

１…流体シミュレーション装置
１０…制御部
１１…シミュレーション受付部
１２…粒子分類部
１３…演算部
１４…出力部
２０…記憶部
２１…粒子情報
２２…逐次結果情報
３０…近傍リスト
４０…自由流体粒子
４１…内部粒子
４２…表面粒子
５０…壁境界粒子
５１…固定流体粒子
５２…圧力境界粒子
５３…非使用粒子
４０１…ＣＰＵ
４０２…入力装置
４０３…モニタ
４０４…出力装置
４０５…媒体読取装置
４０６…インタフェース装置
４０７…通信装置
４０８…ＲＡＭ
４０９…ハードディスク装置
４１０…バス
４１１…プログラム
４１２…各種データ
Ｃ１、Ｃ２…ケース
ｉ、ｉａ、ｉｂ、ｊ、ｊ’、ｊ’’…粒子

Claims (4)

1. 所定の時間における各粒子の粒子データに含まれる位置情報をもとに、壁との境界にかかる壁境界粒子の中から、流体粒子との距離が所定値以下の粒子を第１粒子とし、前記第１粒子との距離が所定値以下の粒子を第２粒子として抽出し、
   抽出した前記第１粒子および前記第２粒子をもとに、各粒子の圧力を求める圧力ポアソン方程式の境界条件を設定して各粒子にかかる圧力を求め、
   求めた前記圧力に基づいて次の時間における各粒子の粒子データを算出する、
   処理をコンピュータに実行させることを特徴とする流体シミュレーションプログラム。
2. 前記求める処理は、前記第１粒子を速度を所定値とする流体粒子とし、前記第２粒子の圧力を０として前記境界条件を設定する、
   請求項１に記載の流体シミュレーションプログラム。
3. 所定の時間における各粒子の粒子データに含まれる位置情報をもとに、壁との境界にかかる壁境界粒子の中から、流体粒子との距離が所定値以下の粒子を第１粒子とし、前記第１粒子との距離が所定値以下の粒子を第２粒子として抽出し、
   抽出した前記第１粒子および前記第２粒子をもとに、各粒子の圧力を求める圧力ポアソン方程式の境界条件を設定して各粒子にかかる圧力を求め、
   求めた前記圧力に基づいて次の時間における各粒子の粒子データを算出する、
   処理をコンピュータが実行することを特徴とする流体シミュレーション方法。
4. 所定の時間における各粒子の粒子データに含まれる位置情報をもとに、壁との境界にかかる壁境界粒子の中から、流体粒子との距離が所定値以下の粒子を第１粒子とし、前記第１粒子との距離が所定値以下の粒子を第２粒子として抽出する分類部と、
   抽出した前記第１粒子および前記第２粒子をもとに、各粒子の圧力を求める圧力ポアソン方程式の境界条件を設定して各粒子にかかる圧力を求め、求めた前記圧力に基づいて次の時間における各粒子の粒子データを算出する算出部と、
   を有することを特徴とする流体シミュレーション装置。

Images