(57)【要約】      （修正有）【課題】粒子法における解の収束性を確保することができる流体シミュレーションプログラム、流体シミュレーション方法および流体シミュレーション装置を提供する。【解決手段】流体シミュレーションプログラムは、抽出する処理と、求める処理と、算出する処理とをコンピュータに実行させる。抽出する処理は、所定の時間における各粒子の粒子データに含まれる位置情報をもとに、壁との境界にかかる壁境界粒子の中から、流体粒子との距離が所定値以下の粒子を第１粒子とし、第１粒子との距離が所定値以下の粒子を第２粒子として抽出する。求める処理は、抽出した第１粒子および第２粒子をもとに、各粒子の圧力を求める圧力ポアソン方程式の境界条件を設定して各粒子にかかる圧力を求める。算出する処理は、求めた圧力に基づいて次の時間における各粒子の粒子データを算出する。【選択図】図４
2019070597.tif 000002
【特許請求の範囲】
【請求項１】
所定の時間における各粒子の粒子データに含まれる位置情報をもとに、壁との境界にかかる壁境界粒子の中から、流体粒子との距離が所定値以下の粒子を第１粒子とし、前記第１粒子との距離が所定値以下の粒子を第２粒子として抽出し、
抽出した前記第１粒子および前記第２粒子をもとに、各粒子の圧力を求める圧力ポアソン方程式の境界条件を設定して各粒子にかかる圧力を求め、
求めた前記圧力に基づいて次の時間における各粒子の粒子データを算出する、
処理をコンピュータに実行させることを特徴とする流体シミュレーションプログラム。
【請求項２】
前記求める処理は、前記第１粒子を速度を所定値とする流体粒子とし、前記第２粒子の圧力を０として前記境界条件を設定する、
請求項１に記載の流体シミュレーションプログラム。
【請求項３】
所定の時間における各粒子の粒子データに含まれる位置情報をもとに、壁との境界にかかる壁境界粒子の中から、流体粒子との距離が所定値以下の粒子を第１粒子とし、前記第１粒子との距離が所定値以下の粒子を第２粒子として抽出し、
抽出した前記第１粒子および前記第２粒子をもとに、各粒子の圧力を求める圧力ポアソン方程式の境界条件を設定して各粒子にかかる圧力を求め、
求めた前記圧力に基づいて次の時間における各粒子の粒子データを算出する、
処理をコンピュータが実行することを特徴とする流体シミュレーション方法。
【請求項４】
所定の時間における各粒子の粒子データに含まれる位置情報をもとに、壁との境界にかかる壁境界粒子の中から、流体粒子との距離が所定値以下の粒子を第１粒子とし、前記第１粒子との距離が所定値以下の粒子を第２粒子として抽出する分類部と、
抽出した前記第１粒子および前記第２粒子をもとに、各粒子の圧力を求める圧力ポアソン方程式の境界条件を設定して各粒子にかかる圧力を求め、求めた前記圧力に基づいて次の時間における各粒子の粒子データを算出する算出部と、
を有することを特徴とする流体シミュレーション装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明の実施形態は、流体シミュレーションプログラム、流体シミュレーション方法および流体シミュレーション装置に関する。
【背景技術】
【０００２】
従来、流体や弾性体などの連続体（以下、流体）の運動を数値的に計算する流体解析手法の一つとして、連続体を粒子として離散化し、粒子分布で表現する粒子法がある。この粒子法としては、Smoothed  Particle  Hydrodynamics（SPH）法や  Moving  Particle  Semi-implicitもしくはMoving  Particle  Simulation（MPS）法が知られている。
【０００３】
これら、粒子分布を用いて解析対象を表現する手法においては、ある粒子から影響半径と呼ばれる距離（ｈ）以内にある粒子を近傍粒子としてその粒子の情報を用いる。一例として、ＳＰＨ法で非圧縮流体の運動方程式として以下の式（１）を解くことを考える。
【０００４】
【数１】
2019070597.tif 000003
【０００５】
ここで、左辺は流体要素の加速度、右辺第一項は圧力勾配項、第二項は外力項である。これをＳＰＨ法を用いて離散化する。ｘｉ，ｖｉ，ρｉ，ｍｉ，ｐｉをそれぞれ、粒子ｉの位置ベクトル、速度ベクトル、密度、質量、圧力とすると、時間発展について以下の式（２ａ）、（２ｂ）のとおりに計算することができる。
【０００６】
【数２】
2019070597.tif 000004
【０００７】
ここで、Δｔはシミュレーションの時間刻み、上付き添え字ｎは時間ステップ数である。また、角括弧＜＞はＳＰＨ法による重ね合わせ評価を表し、次の式（３）のとおりである。
【０００８】
【数３】
2019070597.tif 000005
【０００９】
ここで、和の対象ｊは粒子ｉから半径ｈ以内に存在する粒子(近傍粒子)のインデックスを渡るものとする。また、以後ＳＰＨ法の重ね合わせの和をとるときはすべて同様である。Ｗはカーネル関数で、粒子の分布から連続場を構成するのに用い、以下の式（４）の３次のスプライン関数等がよく使われる。ここで、βは規格化係数であり、解析対象が三次元の場合はπh3  である。
【００１０】
【数４】
2019070597.tif 000006
【００１１】
各粒子の圧力ｐｉは非圧縮条件∇・ｖｉ＝０を満たすように決めるため、式（５）の圧力ポアソン方程式を陰的に解くことで求める。
【００１２】
【数５】
2019070597.tif 000007
【００１３】
ここで、ρ０は物質の基準密度、ｖｉ＊は外力項によって更新された中間状態の速度であり、次の式（６）で表現される。
【００１４】
【数６】
2019070597.tif 000008
【００１５】
式（５）の左辺は圧力のラプラシアンであり、次の式（７）で離散化される。
【００１６】
【数７】
2019070597.tif 000009
【００１７】
また、式（５）の右辺の速度発散は、次の式（８）で離散化される。
【００１８】
【数８】
2019070597.tif 000010
【００１９】
このように離散化することで、式（７）の圧力ポアソン方程式は、次の式（９）、（１０）として表すことができ、次の式（１１）という線形方程式に帰着する。
【００２０】
【数９】
2019070597.tif 000011
【００２１】
【数１０】
2019070597.tif 000012
【００２２】
上記の圧力ポアソン方程式は、ある粒子の近傍に十分な粒子がない時に、その粒子を自由表面とみなし、自由表面の粒子に対して圧力を０と定めて適切な境界条件を与えることで、共役勾配法等の反復法により解くことができる。
【先行技術文献】
【特許文献】
【００２３】
【特許文献１】特開２０１０−２４３２９３号公報
【特許文献２】特開２０１６−１２５９３４号公報
【非特許文献】
【００２４】
【非特許文献１】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）などの半導体メモリ素子、または、ハードディスク、光ディスクなどの記憶装置である。記憶部２０は、粒子情報２１と、逐次結果情報２２とを有する。
【００３４】
粒子情報２１は、流体シミュレーションの設定条件として入力される入力データであり、主には流体シミュレーションにかかる各粒子の情報である。この粒子情報２１における入力データは、例えば計算パラメータと、３種類の各粒子の粒子データとを有する。
【００３５】
計算パラメータは、シミュレーション対象となる流体の物性値やシミュレーション環境などを示すパラメータである。例えば、計算パラメータとしては、各粒子が影響を及ぼす影響半径、各粒子に及ぶ外力ベクトル、シミュレーションの時間刻み幅、シミュレーション対象の物質の基準密度、一粒子の質量、ゼロ除算防止パラメータηなどを含む。
【００３６】
粒子データは、自由流体粒子、壁境界粒子、壁内部粒子の３種類を含む。図２−１は、自由流体粒子の入力データの一例を示す説明図である。図２−２は、壁境界粒子の入力データの一例を示す説明図である。図２−３は、壁内部粒子の入力データの一例を示す説明図である。
【００３７】
自由流体粒子は、シミュレーションの解析対象となる流体を表す粒子である。すなわち、自由流体粒子は、流体粒子の一例である。図２−１に示すように、自由流体粒子にかかる自由流体粒子入力データ１００は、各粒子（ｉ、ｉ＋１、…）の位置座標１０１、速度１０２密度１０３および質量１０４を有する。なお、位置座標１０１および速度１０２については、入力データを初期値として、時間ステップごとに計算された位置および速度の値で更新される。
【００３８】
壁境界粒子は、自由流体粒子の移動が妨げられる固定境界表面に配置し、圧力ポアソン方程式の境界条件を表すための粒子である。図２−２に示すように、壁境界粒子にかかる壁境界粒子入力データ２００は、各粒子（ｉ、ｉ＋１、…）の位置座標２０１、密度２０２および質量２０３を有する。
【００３９】
壁内部粒子は、固定境界内部に配置し、速度発散と圧力勾配を計算する際に用いる粒子である。図２−３に示すように、壁内部粒子にかかる壁内部粒子入力データ３００は、各粒子（ｉ、ｉ＋１、…）の位置座標３０１、法線ベクトル３０２、固定境界面までの距離３０３、密度３０４および質量３０５を有する。なお、これらの粒子の密度、質量については、計算パラメータに含まれる値を使用してもよい。
【００４０】
逐次結果情報２２は、流体シミュレーションにおいて、各粒子の時間ステップごとのシミュレーション結果を示す中間データである。図３−１は、自由流体粒子の中間データの一例を示す説明図である。図３−２は、壁境界粒子の中間データの一例を示す説明図である。図３−３は、壁内部粒子の中間データの一例を示す説明図である。
【００４１】
図３−１に示すように、自由流体粒子にかかる自由流体粒子中間データ１１０は、各粒子（ｉ、ｉ＋１、…）の圧力１１１、中間速度１１２、速度発散１１３および境界フラグ１１４を有する。ここで、境界フラグ１１４は、各自由流体粒子について、内部粒子または表面粒子の分類を示すフラグである。
【００４２】
図３−２に示すように、壁境界粒子にかかる壁境界粒子中間データ２１０は、各粒子（ｉ、ｉ＋１、…）の速度発散２１１、圧力２１２および境界フラグ２１３を有する。ここで、境界フラグ２１３は、各壁境界粒子について、固定流体粒子、圧力境界粒子または非使用粒子の分類を示すフラグである。
【００４３】
図３−３に示すように、壁内部粒子にかかる壁内部粒子中間データ３１０は、各粒子（ｉ、ｉ＋１、…）の仮想マーカ位置３１１、仮想マーカ速度３１２、仮想マーカ加速度３１３、仮想マーカ圧力３１４、中間速度３１５および圧力３１６を有する。ここで、仮想マーカ位置３１１は、壁内部粒子の固定境界面に対する鏡像位置であり、仮想マーカ速度３１２、仮想マーカ加速度３１３および仮想マーカ圧力３１４は、鏡像位置における速度、加速度および圧力の仮想値である。
【００４４】
シミュレーション受付部１１は、流体シミュレーションの実行要求を受け付ける。例えば、シミュレーション受付部１１は、計算パラメータ、自由流体粒子入力データ１００、壁境界粒子入力データ２００および壁内部粒子入力データ３００などの粒子情報２１を含む実行要求を受け付ける。次いで、シミュレーション受付部１１は、受け付けた粒子情報２１を記憶部２０に格納し、シミュレーション受付部１１、粒子分類部１２および出力部１４による流体シミュレーションの処理を開始する。
【００４５】
粒子分類部１２は、計算パラメータの時間刻み幅による時間ステップごとに、粒子情報２１における各粒子の位置情報（位置座標１０１、２０１、３０１）をもとに、各粒子の分類を行う。具体的には、粒子分類部１２は、自由流体粒子について、表面粒子と、内部粒子とに分類する。また、粒子分類部１２は、壁境界粒子について、固定流体粒子と、圧力境界粒子と、非使用粒子とに分類する。
【００４６】
演算部１３は、粒子分類部１２による各粒子の分類結果をもとに、各粒子の圧力を求める圧力ポアソン方程式の境界条件を設定し、各粒子にかかる圧力を求める。具体的には、演算部１３は、壁境界粒子について、固定流体粒子として分類（抽出）された粒子については速度を所定値として固定する流体粒子とし、圧力ポアソン方程式を解く際の自由度に含めるようにする。また、演算部１３は、壁境界粒子について、圧力境界粒子として分類（抽出）された粒子については圧力を０とする境界条件を設定する。次いで、演算部１３は、設定した条件で圧力ポアソン方程式により各粒子にかかる圧力を算出する。次いで、演算部１３は、求めた圧力に基づいて次の時間における各粒子の粒子データ（例えば速度および位置）を算出する。
【００４７】
出力部１４は、演算部１３の算出結果を出力する。具体的には、出力部１４は、時間ステップごとに算出された算出結果を逐次結果情報２２として記憶部２０に格納する。なお、出力部１４は、シミュレーションの開始から終了までの時間ステップごとに算出された逐次結果情報２２に基づくシミュレーション結果を、ディスプレイなどの表示装置への表示出力やファイル出力してもよい。
【００４８】
ここで、流体シミュレーション装置１における流体シミュレーションの処理内容について、詳細を説明する。図４は、実施形態にかかる流体シミュレーション装置１の動作例を示すフローチャートである。
【００４９】
図４に示すように、処理が開始されると、流体シミュレーション装置１のシミュレーション受付部１１は、ユーザからの入力操作などを介して、粒子情報２１にかかる入力データの取得を行う（Ｓ１）。
【００５０】
次いで、流体シミュレーション装置１の粒子分類部１２、演算部１３および出力部１４は、粒子情報２１および逐次結果情報２２を参照し、シミュレーションの終了時間に対応するステップ（時間ステップ）数分のループ処理（Ｓ２〜Ｓ１２）を実行する。
【００５１】
Ｓ２に次いで、粒子分類部１２は、粒子情報２１における位置情報（位置座標１０１、２０１、３０１）に基づいて、粒子ごとの、近傍にある粒子の近傍リストを作成する（Ｓ３）。具体的には、粒子分類部１２は、各粒子において、粒子の位置から計算パラメータにおける影響半径内の距離にある粒子を抽出し、近傍の粒子を識別する粒子番号の近傍リストを作成する。
【００５２】
図５は、近傍リストの一例を示す説明図である。図５に示すように、近傍リスト３０は、粒子（…ｉ−１、ｉ、ｉ＋１、ｉ＋２、…）ごとに、影響半径内にある近傍の粒子を並べたものである。例えば、粒子（ｉ）については、粒子（ｊ１、ｊ２、…）が近傍にある粒子のリストとして挙げられている。また、粒子（ｉ＋１）については、粒子（ｊ’１，ｊ’２、…）が近傍にある粒子のリストとして挙げられている。
【００５３】
次いで、演算部１３は、自由流体粒子の各粒子の速度（ｖｉｎ）に対して、外力ベクトル（ｇ）と時間刻み幅（Δｔ）とによる速度更新（式（６））を行い、中間速度（ｖｉ＊）を算出する（Ｓ４）。
【００５４】
次いで、演算部１３は、自由流体粒子に対して、速度発散を算出する（Ｓ５）。これは、仮想マーカ位置３１１における速度評価、壁内部粒子の速度評価、速度発散の算出、という手順で行う。
【００５５】
具体的には、演算部１３は、仮想マーカ位置３１１を壁内部粒子の固定境界面に対する鏡像位置として算出する。例えば、ある壁内部粒子の仮想マーカ位置３１１は、粒子の位置座標３０１、法線ベクトル３０２および固定境界面までの距離３０３を用いて次の式（１２）で算出される。
【００５６】
【数１１】
2019070597.tif 000013
【００５７】
この位置における速度（仮想マーカ速度３１２）は、次の式（１３）によって評価する。ここで、和をとる範囲は自由流体粒子の近傍粒子とする。
【００５８】
【数１２】
2019070597.tif 000014
【００５９】
演算部１３は、壁内部粒子の中間速度３１５を、仮想マーカ位置３１１の仮想マーカ速度３１２からｖｉ＊＝−ｖｉＶＭとして算出する。これは、固定境界面において非すべり条件を課したことに相当する。
【００６０】
次いで、演算部１３は、各自由流体粒子に対して速度発散を式（８）用いて算出する。ここで、和をとる範囲は、自由流体粒子と壁内部粒子の近傍粒子の範囲とする。
【００６１】
次いで、粒子分類部１２は、自由流体粒子の中から表面粒子とする粒子を検出（抽出）する（Ｓ６）。具体的には、粒子分類部１２は、ある粒子に注目したとき、いずれかの方向に近傍粒子（近傍リスト３０にある、その粒子からの距離が影響半径よりも小さい粒子）が存在しない時を表面粒子とする。また、粒子分類部１２は、全ての方向に近傍粒子が存在するときを内部粒子とする。
【００６２】
図６−１は、粒子の近傍を表す空間領域の分類を説明する説明図である。図６−２は、表面粒子と内部粒子の抽出例を説明する説明図である。
【００６３】
図６−１に示すように、２次元を例にして説明すると、先ず、注目する粒子（ｉ）を中心として影響半径内の空間領域をｘ＋、ｘ−、ｙ＋、ｙ−の４領域に分類する。次いで、粒子分類部１２は、影響半径内の粒子がいずれの空間領域に属するかを判断し、粒子が一つも所属しない領域があった場合、粒子（ｉ）を表面粒子であると抽出する。
【００６４】
図６−２に示すように、ケースＣ１では、ｙ＋の領域に近傍粒子が存在しないため、中心の粒子（ｉａ）については表面粒子として抽出される。一方、ケースＣ２では、全ての領域に近傍粒子が存在するため、粒子（ｉｂ）は内部粒子として抽出される。粒子分類部１２は、この判定結果について、自由流体粒子中間データ１１０の境界フラグ１１４に分類結果として保存する。
【００６５】
次いで、粒子分類部１２は、粒子情報２１における位置座標１０１、２０１、３０１に基づく粒子の分類処理により、壁境界粒子を固定流体粒子、圧力境界粒子、非使用粒子に分類する（Ｓ７）。
【００６６】
図７は、粒子の分類処理を例示するフローチャートである。図７に示すように、分類処理が開始されると、粒子分類部１２は、全ての壁境界粒子の境界フラグ２１３を非使用粒子とする（Ｓ２０）。
【００６７】
次いで、粒子分類部１２は、自由流体粒子ごとのループ処理（ｉ）を実行する（Ｓ２１〜Ｓ２８）。ループ処理（ｉ）が開始されると、粒子分類部１２は、自由流体粒子を一つ選び、粒子ｉとする（Ｓ２２）。
【００６８】
次いで、粒子分類部１２は、壁境界粒子ごとのループ処理（ｊ）を実行する（Ｓ２３〜Ｓ２７）。ループ処理（ｊ）が開始されると、粒子分類部１２は、壁境界粒子を一つ選び粒子ｊとする（Ｓ２４）。
【００６９】
次いで、粒子分類部１２は、粒子ｉの位置座標１０１と、粒子ｊの位置座標２０１とを用いて粒子ｉと粒子ｊの距離（ｒ）を求め、粒子ｉと粒子ｊの距離が影響半径以下であるか否かを判定する（Ｓ２５）。
【００７０】
影響半径以下である場合（Ｓ２５：ＹＥＳ）、粒子分類部１２は、粒子ｊの境界フラグ２１３を固定流体粒子とする（Ｓ２６）。影響半径以下でない場合（Ｓ２５：ＮＯ）、粒子分類部１２は、Ｓ２６をスキップして、境界フラグ２１３をそのままとする。
【００７１】
次いで、粒子分類部１２は、まだ選択されていない壁境界粒子があれば、その粒子を粒子ｊとして、ループ処理（ｊ）を繰り返す。また、粒子分類部１２は、ループ処理（ｊ）の後、まだ選択されていない自由流体粒子があれば、その粒子を粒子ｉとして、ループ処理（ｉ）を繰り返す。
【００７２】
次いで、粒子分類部１２は、固定流体粒子ごとのループ処理（ｊ’）を実行する（Ｓ２９〜Ｓ３７）。ループ処理（ｊ’）が開始されると、粒子分類部１２は、境界フラグ２１３を固定流体粒子とされた壁境界粒子を一つ選び粒子ｊ’とする（Ｓ３０）。
【００７３】
次いで、粒子分類部１２は、壁境界粒子ごとのループ処理（ｊ’’）を実行する（Ｓ３１〜Ｓ３６）。ループ処理（ｊ’’）が開始されると、粒子分類部１２は、壁境界粒子を一つ選び粒子ｊ’’とする（Ｓ３２）。
【００７４】
次いで、粒子分類部１２は、粒子ｊ’’の境界フラグ２１３が非使用粒子であるか否かを判定する（Ｓ３３）。非使用粒子でない場合（Ｓ３３：ＮＯ）、粒子分類部１２は、Ｓ３６へ処理をスキップし、次の粒子ｊ’’を選ぶループ処理を行う。
【００７５】
非使用粒子である場合（Ｓ３３：ＹＥＳ）、粒子分類部１２は、粒子ｊ’、ｊ’’の位置座標２０１を用いて粒子ｊ’、ｊ’’の距離（ｒｊ’ｊ’’）を求め、粒子ｉと粒子ｊの距離が影響半径以下であるか否かを判定する（Ｓ３４）。
【００７６】
影響半径以下である場合（Ｓ３４：ＹＥＳ）、粒子分類部１２は、粒子ｊ’’の境界フラグ２１３を圧力境界粒子とする（Ｓ３５）。影響半径以下でない場合（Ｓ３４：ＮＯ）、粒子分類部１２は、Ｓ３５をスキップして、境界フラグ２１３をそのままとする。
【００７７】
次いで、粒子分類部１２は、まだ選択されていない壁境界粒子があれば、その粒子を粒子ｊ’’として、ループ処理（ｊ’’）を繰り返す。また、粒子分類部１２は、ループ処理（ｊ’’）の後、まだ選択されていない固定流体粒子があれば、その粒子を粒子ｊ’として、ループ処理（ｊ’）を繰り返す。
【００７８】
図８は、粒子の分類を説明する説明図である。図８に示すように自由流体粒子４０については、Ｓ６の処理により、内部粒子４１または表面粒子４２に分類される。また、壁境界粒子５０については、Ｓ７の分類処理により、固定流体粒子５１、圧力境界粒子５２および非使用粒子５３のいずれかに分類される。
【００７９】
図４に戻り、Ｓ７に次ぐ処理を説明する。なお、非使用粒子５３については、このステップが終わるまで使用されないので、以後、無視することとする。ここで改めて自由流体粒子４０、固定流体粒子５１、圧力境界粒子５２などに通し番号を振り、以後、その番号を用いて粒子ｉ、粒子ｊ等として参照する。また、演算部１３は、固定流体粒子５１、圧力境界粒子５２の圧力２１２および速度発散２１１を０とする。
【００８０】
Ｓ７に次いで、演算部１３は、各粒子の圧力を求める圧力ポアソン方程式として、式（１１）に示した線形方程式を解くことで各粒子の圧力を求める（Ｓ８）。
【００８１】
具体的には、演算部１３は、作業データとして係数行列Ａと右辺ベクトルｂを用意する。係数行列Ａは粒子数の２乗個の実数値が格納できれば十分であるが、そのほとんどが０値となるので、非ゼロ要素のみを効率よく格納するためにCompressed  Sparse  Row（CSR）形式などを用いてもよい。右辺ベクトルｂは粒子数と同数の実数値を格納できる配列として確保する。演算部１３は、この係数行列Ａと右辺ベクトルｂに対して行列方程式（式１１）を解く求解処理を行う。
【００８２】
図９は、求解処理を例示するフローチャートである。図９に示すように、処理が開始されると、演算部１３は、係数行列Ａｉｊの値を全て０とし、全要素０で初期化する（Ｓ４０）。
【００８３】
次いで、演算部１３は、自由流体粒子または壁境界粒子を一つ選び粒子ｉとする（Ｓ４１）。次いで、演算部１３は、粒子ｉの境界フラグ１１４または境界フラグ２１３が内部粒子４１または固定流体粒子５１であるか否かを判定する（Ｓ４２）。
【００８４】
内部粒子４１または固定流体粒子５１でない場合（Ｓ４２：ＮＯ）、演算部１３は、係数行列Ａｉｊの要素を１とする（Ｓ４９）。次いで、演算部１３は、右辺ベクトルｂｉの要素の値を粒子ｉの圧力ｐｉ（圧力１１１または圧力２１２）とし（Ｓ５０）、Ｓ５１へ処理を進める。
【００８５】
内部粒子４１または固定流体粒子５１である場合（Ｓ４２：ＹＥＳ）、演算部１３は、粒子ｉの近傍リスト３０に含まれる粒子ごとのループ処理（ｊ）を実行する（Ｓ４３〜Ｓ４７）。ループ処理（ｊ）が開始されると、演算部１３は、粒子ｉの近傍リスト３０に含まれる粒子を一つ選び、粒子ｊとする（Ｓ４４）。
【００８６】
次いで、演算部１３は、係数行列Ａｉｊの要素について、式（９）の右辺上によって値を算出する（Ｓ４５）。ここで、演算部１３は、各粒子の位置座標１０１または位置座標２０１、密度１０３または密度２０２、質量１０４または質量２０３を用いる。次いで、演算部１３は、係数行列Ａｉｉの要素に、Ｓ４５で算出した値を加算する（Ｓ４６）。
【００８７】
次いで、演算部１３は、粒子ｉの近傍リスト３０に含まれる粒子の中でまだ選択されていない粒子があれば、その粒子を粒子ｊとして、ループ処理（ｊ）を繰り返す。
【００８８】
粒子ｉの近傍粒子の全てに対してＳ４３〜Ｓ４７の処理が済んだ後、演算部１３は、右辺ベクトルｂｉの値を式（１０）を用いて算出する（Ｓ４８）。ここでは、各粒子の速度発散１１３または速度発散２１１を用いる。
【００８９】
次いで、演算部１３は、まだ選んでいない内部粒子４１か固定流体粒子５１があるか否かを判定する（Ｓ５１）。まだ選んでいない内部粒子４１か固定流体粒子５１がある場合（Ｓ５１：ＹＥＳ）、演算部１３は、Ｓ４１へ処理を戻し、新たな粒子ｉを選ぶ。
【００９０】
まだ選んでいない内部粒子４１か固定流体粒子５１がない場合（Ｓ５１：ＮＯ）、演算部１３は、ここまでの処理で得られた係数行列Ａｉｊと右辺ベクトルｂおよび各粒子の圧力１１１または圧力２１２をｐとし、行列方程式を解く（Ｓ５２）。具体的には、演算部１３は、行列方程式Ａｉｊｐｊ＝ｂｉを解くことで、圧力を求める。この求解には、例えば共役勾配法などの解法を用いればよい。
【００９１】
図４に戻り、Ｓ８に次いで、演算部１３は、Ｓ８の求解処理による算出結果をもとに、各粒子について圧力勾配による加速度（速度）の算出を行う（Ｓ９）。具体的には、演算部１３は、Ｓ８により各自由流体粒子の圧力ｐｉが求められたところで、次に壁内部粒子の圧力を求める。
【００９２】
演算部１３は、壁内部粒子の圧力を求めるために、先ず仮想マーカ加速度３１３と仮想マーカ圧力３１４を次の式（１４）、（１５）によって求める。
【００９３】
【数１３】
2019070597.tif 000015
【００９４】
【数１４】
2019070597.tif 000016
【００９５】
次に、演算部１３は、式（１４）、（１５）によって求めた値を用いて、壁内部粒子の圧力ｐｉｗａｌｌを次の式（１６）によって算出する。
【００９６】
【数１５】
2019070597.tif 000017
【００９７】
このようにして自由流体粒子と壁内部粒子の圧力が求められたところで、演算部１３は、式（３）を用いて圧力勾配を求め、圧力勾配により粒子の加速度、速度を更新する。なお、式（３）の和については、自由流体粒子と壁内部粒子との近傍をわたるものとする。
【００９８】
次いで、演算部１３は、速度・加速度の時間積分を行い、各粒子の位置を更新して時間を１ステップ進める（Ｓ１０）。具体的には、演算部１３は、求めた速度に１ステップ分の時間をかけ合わせて各粒子の移動距離を求め、各粒子の位置（位置座標１０１）を更新することで１ステップ分の計算が完了する。
【００９９】
次いで、出力部１４は、計算結果を逐次結果情報２２としてファイル出力する（Ｓ１１）。流体シミュレーション装置１は、このＳ２〜Ｓ１２のループ処理を繰り返すことで、シミュレーションの開始から終了までの、計算パラメータの時間刻み幅による時間ステップごとのシミュレーション結果を得る。
【０１００】
以上のように、流体シミュレーション装置１は、粒子分類部１２と、演算部１３とを有する。粒子分類部１２は、各粒子の時間ステップごとの粒子データに含まれる位置情報をもとに、壁境界粒子５０の中から、自由流体粒子４０との距離が所定値以下の粒子を固定流体粒子５１として抽出する。また、粒子分類部１２は、固定流体粒子５１との距離が所定値以下の粒子を圧力境界粒子５２として抽出する。演算部１３は、抽出した固定流体粒子５１および圧力境界粒子５２をもとに、各粒子の圧力を求める圧力ポアソン方程式の境界条件を設定して各粒子にかかる圧力を求める。次いで、演算部１３は、求めた圧力に基づいて次の時間ステップにおける各粒子の粒子データを算出する。このように壁境界粒子５０より抽出した固定流体粒子５１および圧力境界粒子５２をもとに各粒子の圧力を求める圧力ポアソン方程式の境界条件を設定することで、流体シミュレーション装置１は、例えば細い隙間に流体が入り込む場合であっても適切な境界条件を設定することができる。このため、流体シミュレーション装置１は、粒子法における解の収束性を確保することができる。
【０１０１】
ここで、適切な境界条件を設定し、境界条件を改善することで、粒子法における解の収束性を確保されることを、解が発散するケースと比較して説明する。図１０−１は、解が発散するケースを説明する説明図である。図１０−２は、境界条件の改善例を説明する説明図である。
【０１０２】
図１０−１に示すように、壁境界粒子５０が一律に設定されるケースでは、狭い隙間に入り込む自由流体粒子４０について、粒子４０ａを取り囲む粒子が多く、内部粒子として判定される。このケースでは、粒子４０ａの近傍にある壁境界粒子５０において圧力ポアソン方程式を解く際の境界条件に含める粒子が存在しないことから、境界条件が設定されずに圧力ポアソン方程式を解くことで解が発散する場合がある。
【０１０３】
これに対し、本実施形態では、壁境界粒子５０の中から、自由流体粒子４０との距離が所定値以下の粒子を固定流体粒子５１として抽出し、固定流体粒子５１との距離が所定値以下の粒子を圧力境界粒子５２として抽出する。そして、壁境界粒子５０より抽出した固定流体粒子５１および圧力境界粒子５２をもとに各粒子の圧力を求める圧力ポアソン方程式の境界条件を設定する。
【０１０４】
例えば、流体シミュレーション装置１は、固定流体粒子５１を速度固定の流体粒子とし、圧力境界粒子５２を圧力を０として圧力ポアソン方程式の境界条件を設定する。このように、流体シミュレーション装置１は、流体の近傍となった固定流体粒子５１を圧力ポアソン方程式を解く際の自由度に含め、固定流体粒子５１近傍の圧力境界粒子５２は圧力０の表面粒子として圧力ポアソン方程式の境界条件とする。これにより、流体シミュレーション装置１は、例えば細い隙間に流体が入る事象を破綻なく計算できる。
【０１０５】
具体的には、図１０−２に示すように、自由流体粒子４０の近傍にある壁境界粒子５０については、例えば速度固定の固定流体粒子５１として扱い、圧力ポアソン方程式の自由度に含める。また、流体の近傍ではないが、固定流体粒子５１の近傍にある壁境界粒子を圧力境界粒子５２（表面粒子）として、圧力０と定めることで、圧力ポアソン方程式の境界条件とする。なお、流体の近傍でも固定流体粒子５１の近傍でもない壁境界粒子５０については、非使用粒子５３として圧力ポアソン方程式の算出時に使用しないものとする。このように、固定流体粒子５１および圧力境界粒子５２をもとに各粒子の圧力を求める圧力ポアソン方程式の境界条件を設定することで、狭い隙間に入り込む自由流体粒子４０がある場合でも、解の収束性を確保することができる。
【０１０６】
なお、図示した各装置の各構成要素は、必ずしも物理的に図示の如く構成されていることを要しない。すなわち、各装置の分散・統合の具体的形態は図示のものに限られず、その全部または一部を、各種の負荷や使用状況などに応じて、任意の単位で機能的または物理的に分散・統合して構成することができる。
【０１０７】
また、流体シミュレーション装置１で行われる各種処理機能は、ＣＰＵ（またはＭＰＵ、ＭＣＵ（Micro  Controller  Unit）等のマイクロ・コンピュータ）上で、その全部または任意の一部を実行するようにしてもよい。また、各種処理機能は、ＣＰＵ（またはＭＰＵ、ＭＣＵ等のマイクロ・コンピュータ）で解析実行されるプログラム上、またはワイヤードロジックによるハードウエア上で、その全部または任意の一部を実行するようにしてもよいことは言うまでもない。また、流体シミュレーション装置１で行われる各種処理機能は、クラウドコンピューティングにより、複数のコンピュータが協働して実行してもよい。
【０１０８】
ところで、上記の実施形態で説明した各種の処理は、予め用意されたプログラムをコンピュータで実行することで実現できる。そこで、以下では、上記の実施形態と同様の機能を有するプログラムを実行するコンピュータ（ハードウエア）の一例を説明する。図１１は、実施形態にかかる流体シミュレーション装置１のハードウエア構成例を示すブロック図である。
【０１０９】
図１１に示すように、流体シミュレーション装置１は、各種演算処理を実行するＣＰＵ４０１と、データ入力を受け付ける入力装置４０２と、モニタ４０３と、音声等の出力装置４０４とを有する。また、流体シミュレーション装置１は、記憶媒体からプログラム等を読み取る媒体読取装置４０５と、各種装置と接続するためのインタフェース装置４０６と、有線または無線により外部機器と通信接続するための通信装置４０７とを有する。また、流体シミュレーション装置１は、各種情報を一時記憶するＲＡＭ４０８と、ハードディスク装置４０９とを有する。また、流体シミュレーション装置１内の各部（４０１〜４０９）は、バス４１０に接続される。
【０１１０】
ハードディスク装置４０９には、上記の実施形態で説明したシミュレーション受付部１１、粒子分類部１２、演算部１３および出力部１４などで各種の処理を実行するためのプログラム４１１が記憶される。また、ハードディスク装置４０９には、プログラム４１１が参照する各種データ４１２（例えば粒子情報２１、逐次結果情報２２等）が記憶される。入力装置４０２は、例えば、流体シミュレーション装置１の操作者から操作情報の入力を受け付ける。モニタ４０３は、例えば、操作者が操作する各種画面を表示する。インタフェース装置４０６は、例えば印刷装置等が接続される。通信装置４０７は、ＬＡＮ（Local  Area  Network）等の通信ネットワークと接続され、通信ネットワークを介した外部機器との間で各種情報をやりとりする。
【０１１１】
ＣＰＵ４０１は、ハードディスク装置４０９に記憶されたプログラム４１１を読み出して、ＲＡＭ４０８に展開して実行することで、シミュレーション受付部１１、粒子分類部１２、演算部１３および出力部１４にかかる各種の処理を行う。なお、プログラム４１１は、ハードディスク装置４０９に記憶されていなくてもよい。例えば、流体シミュレーション装置１が読み取り可能な記憶媒体に記憶されたプログラム４１１を読み出して実行するようにしてもよい。流体シミュレーション装置１が読み取り可能な記憶媒体は、例えば、ＣＤ−ＲＯＭやＤＶＤディスク、ＵＳＢ（Universal  Serial  Bus）メモリ等の可搬型記録媒体、フラッシュメモリ等の半導体メモリ、ハードディスクドライブ等が対応する。また、公衆回線、インターネット、ＬＡＮ等に接続された装置にこのプログラム４１１を記憶させておき、流体シミュレーション装置１がこれらからプログラム４１１を読み出して実行するようにしてもよい。
【０１１２】
以上の実施形態に関し、さらに以下の付記を開示する。
【０１１３】
（付記１）所定の時間における各粒子の粒子データに含まれる位置情報をもとに、壁との境界にかかる壁境界粒子の中から、流体粒子との距離が所定値以下の粒子を第１粒子とし、前記第１粒子との距離が所定値以下の粒子を第２粒子として抽出し、
抽出した前記第１粒子および前記第２粒子をもとに、各粒子の圧力を求める圧力ポアソン方程式の境界条件を設定して各粒子にかかる圧力を求め、
求めた前記圧力に基づいて次の時間における各粒子の粒子データを算出する、
処理をコンピュータに実行させることを特徴とする流体シミュレーションプログラム。
【０１１４】
（付記２）前記求める処理は、前記第１粒子を速度を所定値とする流体粒子とし、前記第２粒子の圧力を０として前記境界条件を設定する、
付記１に記載の流体シミュレーションプログラム。
【０１１５】
（付記３）所定の時間における各粒子の粒子データに含まれる位置情報をもとに、壁との境界にかかる壁境界粒子の中から、流体粒子との距離が所定値以下の粒子を第１粒子とし、前記第１粒子との距離が所定値以下の粒子を第２粒子として抽出し、
抽出した前記第１粒子および前記第２粒子をもとに、各粒子の圧力を求める圧力ポアソン方程式の境界条件を設定して各粒子にかかる圧力を求め、
求めた前記圧力に基づいて次の時間における各粒子の粒子データを算出する、
処理をコンピュータが実行することを特徴とする流体シミュレーション方法。
【０１１６】
（付記４）前記求める処理は、前記第１粒子を速度を所定値とする流体粒子とし、前記第２粒子の圧力を０として前記境界条件を設定する、
付記３に記載の流体シミュレーション方法。
【０１１７】
（付記５）所定の時間における各粒子の粒子データに含まれる位置情報をもとに、壁との境界にかかる壁境界粒子の中から、流体粒子との距離が所定値以下の粒子を第１粒子とし、前記第１粒子との距離が所定値以下の粒子を第２粒子として抽出する分類部と、
抽出した前記第１粒子および前記第２粒子をもとに、各粒子の圧力を求める圧力ポアソン方程式の境界条件を設定して各粒子にかかる圧力を求め、求めた前記圧力に基づいて次の時間における各粒子の粒子データを算出する算出部と、
を有することを特徴とする流体シミュレーション装置。
【０１１８】
（付記６）前記算出部は、前記第１粒子を速度を所定値とする流体粒子とし、前記第２粒子の圧力を０として前記境界条件を設定する、
付記５に記載の流体シミュレーション装置。
【符号の説明】
【０１１９】
１…流体シミュレーション装置
１０…制御部
１１…シミュレーション受付部
１２…粒子分類部
１３…演算部
１４…出力部
２０…記憶部
２１…粒子情報
２２…逐次結果情報
３０…近傍リスト
４０…自由流体粒子
４１…内部粒子
４２…表面粒子
５０…壁境界粒子
５１…固定流体粒子
５２…圧力境界粒子
５３…非使用粒子
４０１…ＣＰＵ
４０２…入力装置
４０３…モニタ
４０４…出力装置
４０５…媒体読取装置
４０６…インタフェース装置
４０７…通信装置
４０８…ＲＡＭ
４０９…ハードディスク装置
４１０…バス
４１１…プログラム
４１２…各種データ
Ｃ１、Ｃ２…ケース
ｉ、ｉａ、ｉｂ、ｊ、ｊ’、ｊ’’…粒子
【図１】
2019070597.tif 000018
【図２−１】
2019070597.tif 000019
【図２−２】
2019070597.tif 000020
【図２−３】
2019070597.tif 000021
【図３−１】
2019070597.tif 000022
【図３−２】
2019070597.tif 000023
【図３−３】
2019070597.tif 000024
【図４】
2019070597.tif 000025
【図５】
2019070597.tif 000026
【図６−１】
2019070597.tif 000027
【図６−２】
2019070597.tif 000028
【図７】
2019070597.tif 000029
【図８】
2019070597.tif 000030
【図９】
2019070597.tif 000031
【図１０−１】
2019070597.tif 000032
【図１０−２】
2019070597.tif 000033
【図１１】
2019070597.tif 000034