(57)【要約】【課題】解析空間中に置かれた移動物体周りの流れ解析を、境界埋め込み法にて解析する場合に、各計算セルの体積と各計算セルを占める物体体積の比である物体体積分率と、移動物体が流れに与える影響を反映させるため相互作用力を計算する必要がある。これらの計算の精度向上と負荷低減の両立が困難であった。【解決手段】計算セル４１００を、物体表面４２００との関係から、計算セル全体が物体の外部に位置する計算セル（以下、外部セル）４１０７，計算セル全体が物体内部に含まれる計算セル（以下、内部セル）４１０８，計算セル内部に物体境界を含む計算セル（境界面上セル）４１０６の３種類に分類する。境界面上セル４１０６に分類される計算セルにおける物体体積分率及び相互作用力を、物体内部に散らばった仮想粒子の体積，移動速度に基づき計算する。【選択図】  図２０
2011070333.tif 000002
【特許請求の範囲】
【請求項１】
解析空間に複数散らばった格子点の位置座標により表される複数の格子点データと、
前記複数の格子点データを用いて表される複数のセル表面データと、
前記複数のセル表面データにより表される計算セルデータとから構成される計算格子データと、
解析空間に置かれた物体を表すために物体表面に複数散らばった点の位置座標により表される物体点データと、
前記複数の物体点データにより表される物体微小面データから構成される物体データを記憶する手段とを備える、
前記物体周りの流れを解析するための流体解析装置において、
前記解析空間に置かれた物体の内部に分布した複数の仮想粒子の位置座標と物体体積及び格子データに基づき決定した体積データで構成される仮想粒子データを記憶する手段を備え、
前記物体点データ及び前記物体微小面データに基づき前記格子点が移動物体の内外のいずれにあるか判定する格子点判定手段を備え、
各計算セル中に含まれる物体体積と前記各計算セルの体積の比である物体の体積分率を、前記格子点判定手段で物体内部にあると判定した格子点である内部点と物体外部にあると判定した格子点である外部点の両方を用いて構成される計算セルである境界面上セルにおいては、前記境界面上セルに所属する仮想粒子の体積データ基づき計算すること、
を特徴とする流体解析装置。
【請求項２】
請求項１の流体解析装置で、
境界面上セルに分類される一つの計算セルの体積データＶc、前記境界面上セルに分類された一つの計算セルに所属する仮想粒子数Ｎpc、前記仮想粒子の体積Ｖpに基づき、前記境界面上セルの物体体積分率θcを求めること、
を特徴とする流体解析装置。
【請求項３】
請求項１及び請求項２の流体解析装置で、
前記内部点でのみ構成される計算セルである内部セルの体積の総和ＶsumＣを求める計算手段を備え、
前記境界面上セルに分類される計算セルに所属する仮想粒子の総和Ｎtを求める計算手段を備え、
前記物体の体積Ｖmと前記ＶsumＣと前記Ｎtに基づき、前記境界面上格子に分類される計算セルに所属する各仮想粒子の体積Ｖpを式（Ａ）で求めること
Ｖp＝ＶsumＣ／Ｎｔ                        式（Ａ）
を特徴とする流体解析装置。
【請求項４】
請求項３の流体解析装置で、
前記境界面上セルに分類される一つの計算セルに所属する仮想粒子数Ｎpc、前記境界面上セルに分類された一つの計算セルの体積データＶc、前記各仮想粒子の体積Ｖpに基づき前記境界面上セルの物体体積分率θcを式（Ｂ）で求めること
θc＝Ｎpc＊Ｖp／Ｖc                        式（Ｂ）
を特徴とする流体解析装置。
【請求項５】
請求項１から請求項４の流体解析装置で、
前記内部点のみで構成される計算セルである内部セルに分類される各計算セルにおける前記体積分率を１と計算し、かつ、前記外部セルに分類される各計算セルにおいて物体の体積分率を０とすること、
を特徴とする流体解析装置。
【請求項６】
請求項１から請求項５の流体解析装置で、
前記物体が、解析空間中を任意の速度で移動すること、
を特徴とする流体解析装置。
【請求項７】
請求項１から請求項６の流体解析装置で、
前記物体内部に発生させる前記仮想粒子の個数Ｍpを入力する手段と記憶する手段を備え、
前記仮想粒子を発生させる範囲を示す物体表面からの距離Ｌpを入力する手段と記憶する手段を備え、
前記Ｍpと前記Ｌpに基づき仮想粒子を発生すること
を特徴とする流体解析装置。
【請求項８】
請求項１から請求項７の流体解析装置で、
前記物体表面と、前記物体内部に発生させる前記仮想粒子を表示する表示手段を備えること、
を特徴とする流体解析装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、パーソナルコンピュータや大型計算機を用いた流体の数値シミュレーション解析に係り、特に、解析空間に置かれた任意の境界形状を持つ移動物体周りの流れの数値シミュレーションを実行する流体解析装置に関する。
【背景技術】
【０００２】
コンピュータ及び数値解析手法の発達に伴い、流れや燃焼の数値シミュレーションを用いたガソリンエンジン等のレシプロエンジン（以下、エンジン）の筒内現象把握、及びエンジン開発が実施されている。
【０００３】
エンジンには、外部から吸入する空気量を制御するための吸気弁，外部へ排出するガス量を制御するための排気弁，仕事を取り出すためのピストンが設けられており、エンジンの運転中は、これらが移動，停止を繰り返す。このため、エンジン内部の空間形状は時々刻々と変化する。エンジン内部を数値シミュレーションにて解析する場合には、ピストンやバルブの移動に伴う空間形状の変化を何らかの方法で再現する必要がある。
【０００４】
このような空間形状の変化を再現する手法の一つとして境界条件埋め込み法（Immersed Boundary Method、以下、ＩＢ法）が考案されている。流体の数値シミュレーションでは、解析空間を複数の多面体（以下、計算セル）に分割し解析を行う。ＩＢ法は、計算セルとは別に移動物体（例えば、吸気弁，排気弁，ピストン）の表面データを用い、解析空間中でこの表面データを移動させることで、解析空間の形状変化を表現する手法である。移動物体による流れの変化は、流体の運動方程式に相互作用力項を加えることで反映させる。具体的には、式（１）を解くことで流れ解析を行う。
【０００５】
【数１】
2011070333.tif 000003
【０００６】
ここで、ρは流体の密度、ｐは圧力、τはせん断応力テンソル、ｖeｃ(ｕ)は平均速度ベクトル、ｖeｃ(Ｆs)が相互作用力ベクトルである。平均速度ベクトルｖeｃ(ｕ)、相互作用力項ｖeｃ(Ｆs)は式（２），式（３）で評価される。
【０００７】
【数２】
2011070333.tif 000004
【０００８】
【数３】
2011070333.tif 000005
【０００９】
ここで、θは計算セルを占める物体の体積分率（以下、物体体積分率）、ｖeｃ(ｕF)は流体の速度ベクトル、ｖeｃ(ｕP)は物体の速度ベクトル、Δｔはシミュレーションの条件として設定する時間刻み幅であり、解析の安定性や求める計算精度を実現するように決定する。物体の体積分率θは、各計算セルの体積をＶkとし、該計算セルを占める物体の体積をＶMとすると、式（４）で定義される。
【００１０】
【数４】
2011070333.tif 000006
【００１１】
式（４）で定義される物体体積分率は、与えられた移動物体に関するデータ（以下、移動体データ）と計算セルに関するデータ（以下、計算格子データ）に基づき計算される。
【００１２】
ＩＢ法は、計算セル作成の工数を削減できる点がメリットとなる反面、計算格子データと移動体データに基づく体積分率計算の負荷が大きい点がデメリットとなる。このデメリットを解決する手段として、物体形状に基づき配置された仮想粒子の情報に基づき相互作用力項を計算する手法がある（例えば、非特許文献１参照）。
【先行技術文献】
【非特許文献】
【００１３】
【非特許文献１】S.V.Apte、“A multilevel formulation to simulate particulate flows”、CTR Annual Research Report 2004、pp201-208(2004)
【発明の概要】
【発明が解決しようとする課題】
【００１４】
文献によれば、移動物体の内部に複数の仮想粒子を分布させ、移動物体の形状を近似、近似に用いた仮想粒子データに基づき、各計算セルにおける物体の体積分率を計算する方法が報告されている。本手法によれば、幾何的な計算に基づき体積計算をする必要が無くなり、計算負荷低減が図れる。しかし、該手法では、移動物体の境界面の近くにおいて、移動物体が存在しない計算セルで体積分率が０とならない、移動物体で完全に埋め尽くされる格子において、物体の体積分率が１とならない場合がある。これは、形状近似に用いた仮想粒子データのみに基づき計算をしたため、物体の境界の正確な位置情報が失われるために生じると考えられる。この結果、流体のシミュレーション誤差を大きくすることとなる。
【００１５】
本発明は、このような従来技術の課題を解決するものであり、その目的は、移動物体の表面付近に位置する計算セルにおける物体の体積分率を高精度に予測し、これにより、解析時間の精度低下と計算時間の長期化を防ぐ流体解析精度の向上を実現する流体解析装置を提供することである。
【課題を解決するための手段】
【００１６】
上記目的を達成するため、本発明では、解析空間に複数散らばった格子点の位置座標により表される格子点データ、複数の格子点データを用いて表される複数のセル表面データ及び複数のセル表面データにより表される計算セルデータから構成される計算格子データと、解析空間に置かれた物体を表すために物体表面に設けられた複数の物体点データ、複数の物体点データにより表される物体微小面データから構成される物体データを記憶する手段を備える。
【００１７】
前記物体周りの流れを解析するための流体解析装置において、
前記解析空間に置かれた物体の内部に分布した複数仮想粒子の位置座標と物体体積及び格子データに基づき決めた体積データで構成される仮想粒子データを記憶する手段を備え、
前記物体点データ及び前記物体微小面データに基づき前記格子点が移動物体の内外のいずれにあるか判定する格子点判定手段を備え、
各計算セル中に含まれる物体体積と前記各計算セルの体積の比である物体の体積分率を、前記境界面上セルにおいては、前記格子点判定手段で物体内部にあると判定した格子点である内部点と物体外部にあると判定した格子点である外部点の両方で構成される計算セルである境界面上セルに所属する仮想粒子の体積データ基づき計算すること、
を特徴とする流体解析装置
を作成する。
【発明の効果】
【００１８】
本発明によれば、解析空間に設けた複数の計算セルを、内部が移動物体で１００％占められている計算セル（以下、内部セル），内部に移動物体の表面が存在する計算セル（以下、境界面上セル），セル内部に移動物体を含まない計算セル（以下、外部セル）の３種類に分類し、内部セルの物体体積分率を１、外部セルの物体体積分率を０と設定することができる。さらに、境界面上セルの物体体積分率は、移動物体の内部かつ移動物体の表面近くに分布させた仮想粒子の位置及び体積データに基づき計算することができるため、物体体積分率の計算精度を向上できる。この結果、流体解析の精度低下の防止と計算時間長期化の防止を両立できる。
【図面の簡単な説明】
【００１９】
【図１】本発明の実施形態における流体解析装置の全体概略構成図である。
【図２】流体解析装置にて実行する処理の機能ブロック図である。
【図３】流体解析装置にて実行する処理のフローチャートである。
【図４】流れ・物体移動解析手段における処理のフローチャートである。
【図５】移動物体計算手段にて実行する処理の機能ブロック図である。
【図６】移動物体計算手段における処理のフローチャートである。
【図７】格子点区分手段における処理のフローチャートである。
【図８】計算セル区分手段における処理のフローチャートである。
【図９】仮想粒子区分手段における処理のフローチャートである。
【図１０】体積分率・相互作用力計算手段における処理のフローチャートである。
【図１１】物体微小面を用いた移動物体の表現を一部取り出した図である。
【図１２】移動物体データの構成を示した図である。
【図１３】立方体で作成した計算格子の一つを取り出した図である。
【図１４】計算格子データの構成を示した図である。
【図１５】２次元解析空間に分布された仮想粒子を示す図である。
【図１６】仮想粒子データの構成を示した図である。
【図１７】格子点と計算セルの分類を示す図である。
【図１８】仮想粒子発生処理方法を説明する図である。
【図１９】仮想粒子発生条件の入力画面の図である。
【図２０】本発明の実施形態において解析した結果の表示画面である。
【発明を実施するための形態】
【００２０】
図１から図２０を用いて本発明の実施例を説明する。流体解析装置の構成，使用するデータ，処理手段及び処理方法，データ入出力の順に説明する。
【００２１】
図１は、流体解析装置１０００の構成図である。流体解析装置１０００は、キーボード，マウスなどの入力装置１００１，流体解析に用いる格子データ，初期条件，境界条件，移動物体データ，シミュレーション結果などを記録したハードディスクやコンパクトディスク等の外部記憶媒体１００２、シミュレーション結果やシミュレーションを実行するために必要な設定を表示するディスプレイ等の表示装置１００３、および、シミュレーションを実行する解析装置１００４から構成される。解析装置１００４は、入力装置１００１から入力されるデータ、外部記憶媒体１００２との間で入出力されるデータ及び中央処理装置１００６で実行するプログラム等を記憶するための主記憶装置１００５とシミュレーションの処理を記述したプログラムを実行する中央処理装置１００６で構成される。以上を以って流体解析装置を構成する。
【００２２】
続いて、流体解析装置にて使用するデータについて説明する。データは、移動物体を表す移動物体データ、解析空間を分割し生成する計算格子データ、物体が各計算セルを占める割合である物体体積分率の計算に用いる仮想粒子データの順に説明する。
【００２３】
図１１と図１２を用いて移動物体データについて説明する。移動物体データは物体表面の位置と移動速度を表し、かつ、複雑な形状を持つ移動物体に対応する必要がある。このため、移動物体は複数の多角形の集合で表す。以下では物体表面を表すために用いる多角形を物体微小面４２０１と呼ぶ。図１１は移動物体の表面４２００の一部を表示した図である。物体表面４２００は、多数の物体微小面４２０１で構成される。さらに、物体微小面４２０１は、３点の物体点４２０２により構成される。流体解析装置１０００では、これらが数値データの組み合わせで表現される。図１２を用いて、物体点４２０２を表す物体点データ４２０３と物体微小面４２０１を表す物体微小面データ４２０４を説明する。物体点データ４２０３は、各物体点を区別するための物体点番号ＮＢg、物体点の座標，物体点の移動速度で構成される。本実施例では、物体点の総数をＮＢgmと表す。座標は解析空間における位置を示すデータ、物体点の移動速度は、解析空間中を各物体点が移動する速度を示すデータである。物体微小面データ４２０４は、各物体微小面４２０１を区別するための物体面番号ＮＢsと物体微小面４２０１を構成する物体点４２０２の物体点番号で構成される。本実施例では、物体微小面の総数をＮＢsmと表す。物体面番号は整数で与えられるデータである。なお、図１１で三角形を用いて説明したが、微小表面の形状はこれに限るものではなく、四角形や、八角形等の多角形にすることもできる。また、移動物体データの形式は、図１２に示す形式に限るものではない。
【００２４】
続いて、図１３と図１４を用いて計算格子データを説明する。前述の通り、シミュレーションで流体解析を行う際には空間を複数の多面体に分割する。この多面体を計算セルと呼ぶ。図１３には、一例として立方体で作成した計算セル４１００を示している。計算セル４１００は、６面のセル表面４１０１から構成され、各セル表面４１０１は４点の格子点４１０２で構成される。流体解析装置１０００では、これらが数値データの組み合わせで表現される。図１４を用いて、格子点４１０２を表す格子点データ４１０３，セル表面４１０１を表すセル表面データ４１０４，計算セル４１００を表すセルデータ４１０５を説明する。格子点データ４１０３は、各格子点を区別するための格子点番号Ｎg、格子点の座標，物体内部フラグで構成される。本実施例では、格子点の全数をＮgmとする。格子点の座標は、解析空間中における各格子点の位置を示すデータである。図１７を用いて物体内部フラグを説明する。図１７は、２次元の解析空間における計算セル４１００と移動物体の表面４２００が示されている。移動物体の表面４２００との関係から、格子点４１０２は、移動物体の外部にある格子点（以下、外部点）４１１２，移動物体の内部にある格子点（以下、内部点）４１１１の２種類に分類できる。物体内部フラグは、格子点が内部点であるか、外部点であるかを示すデータであり、本実施例では、格子点が内部点であればＴrue、外部点であればＦalseと設定される。セル表面４１０１を表すセル表面データ４１０４は、各セル表面を区別するためのセル表面番号Ｎsと、セル表面を構成する格子点番号Ｎgで構成される。本実施例では、セル表面の総数をＮsmと表す。計算セル４１００を表す計算セルデータ４１０５は、各計算セルを区別するためのセル番号Ｎcと、計算セルを構成するセル表面番号，セルの体積，セル位置フラグ及び各セルを占める物体の体積分率で構成される。本実施例では、セルの総数をＮcmと表す。図１７を用いて、セル位置フラグについて説明する。図１７は、前述の通り２次元の解析空間における計算セル４１００と移動物体の表面４２００が示されている。計算セル４１００は、物体表面４２００との関係から、計算セル全体が物体の外部に位置する計算セル（以下、外部セル）４１０７，計算セル全体が物体内部に含まれる計算セル（以下、内部セル）４１０８，計算セル内部に物体境界を含む計算セル（境界面上セル）４１０６の３種類に分類できる。セル位置フラグは、計算セルが３種類のどれに分類されるかを示すデータであり、本実施例では、計算セルが内部セルであれば０、外部セルであれば１、境界面上セルであれば２と設定される。なお、図１３で立方体セルを用いて説明したが、計算セルの形状はこれに限るものではなく、４面体や、８面体等の多面体にすることもできる。また、計算格子データの形式は、図１４に示す形式に限るものではない。
【００２５】
次に、図１５を用いて体積計算に用いる仮想粒子データについて説明する。図１５には、２次元の解析空間における計算セル４１００，物体表面４２００，仮想粒子４３００が示してある。図１５に示すように、仮想粒子は、移動物体の表面で形成される閉曲面（解析空間が２次元の空間であれば閉曲線）の内側に分布させることが特徴である。仮想粒子は、流体解析装置１０００では、数値データの組み合わせで表現される。仮想粒子４３００を表す仮想粒子データ４３０１は、各仮想粒子を区別するための仮想粒子番号Ｎp，仮想粒子の位置座標，仮想粒子の移動速度，仮想粒子の体積，仮想粒子の所属微小面番号，所属セル番号，物体微小面フラグから構成される。本実施例では、仮想粒子の総数をＮpmと表す。仮想粒子の位置座標及び移動速度は、解析空間に設けた座標軸で与える。所属微小面番号には、物体微小面番号ＮＢsが保存される。例えば、最も距離が近い物体微小面４２０１の番号と設定できる。所属セル番号は、仮想粒子が属する計算セル番号Ｎcが設定される。物体微小面フラグは、仮想粒子が、境界面上セル４１０６に所属するかを示す。境界面上セル４１０６に所属する場合には、Ｔrueが設定され、境界面上セル４１０６以外に所属する場合にはＦalseが設定される。
【００２６】
本実施例では、移動物体を剛体とする。移動物体の移動速度をｖeｃ(ｕP)とし、各物体点、各仮想粒子の移動速度もｖeｃ(ｕP)であるとする。
【００２７】
以上で、データの説明を終わり、続いて処理の説明を行う。処理は、図２と図３により処理の全体像を説明し、その後、各処理の詳細をフローチャートやブロック図を用いて説明する。
【００２８】
図２は流体解析装置１０００で実行する解析機能を表す機能ブロック図である。解析機能は、入力装置１００１や外部記憶媒体１００２から入力されるデータを記憶するための入力データ記憶手段１００，各計算セルにおける物体の体積分率や流体解析を実行する解析手段２００，解析結果を表示装置１００３に表示するための結果表示手段３００で構成される。
【００２９】
入力データ記憶手段１００は、流体解析装置１０００にて解析する物理量（速度，温度，化学種濃度，圧力，乱流等）の初期条件やシミュレーションに必要な設定（例えば、計算時間や時間刻み幅等）、前記物理量の方程式に対する境界条件を記憶する流体解析条件記憶手段１０１，流体解析の解析空間を複数の多面体に分割することで生成される計算格子データを記憶するための計算格子データ記憶手段１０２，移動物体データを記憶するための移動物体データ記憶手段１０３及び仮想粒子データを記憶するための仮想粒子データ記憶手段１０４から構成される。また、解析手段２００は、物体の移動の計算や流体解析条件や計算格子データ，物体体積分率等に基づく流れ場の解析を行う流れ場解析手段２０１，計算格子データ，物体面データ，仮想粒子データに基づき物体体積分率や相互作用力を計算するための移動物体計算手段２０２，仮想粒子データに基づき仮想粒子発生をさせる仮想粒子発生手段２０３から構成される。
【００３０】
図３は、図２に示した解析機能による処理を表すフローチャートである。最初に、ステップＳ１０００にて、解析に必要となるデータを入力装置１００１及び外部記憶媒体１００２から入力する。続いてステップＳ２０００に進み、物体体積分率の計算，流体解析を実行する。続いて、ステップＳ３０００に進み、解析結果をディスプレイ等の表示装置１００３に表示する処理を実行する。
【００３１】
続いて、図３に示したステップＳ２０００の処理を、図４を用いて説明する。図４は、流れ・物体移動解析手段Ｓ２０００の処理をフローチャートで示したものである。ステップＳ２１００にて、解析を開始。ステップＳ２２００にて体積計算に用いる仮想粒子が解析空間に配置されているかを判定する。ステップＳ２２００にて、仮想粒子を配置済み（ＹＥＳ）と判定した場合は、ステップＳ２４００に進む。ステップＳ２２００にて、仮想粒子を未配置（ＮＯ）と判定した場合は、ステップＳ２３００に進み、移動物体内部かつ移動物体の表面近傍に仮想粒子を発生させた後、ステップＳ２４００に進む。ステップＳ２４００では、初期条件に基づき流体解析を開始する。次にステップＳ２５００に進み、時刻ｔにおける物体の体積分率・流れ場の計算を開始する。続いて、ステップＳ２６００に進み、各計算セルを占める移動物体の体積分率や移動物体の影響を反映させる相互作用力を計算し、ステップＳ２７００に進む。ステップＳ２７００では、ステップＳ２６００で計算した物体の体積分率や相互作用力及び入力した計算格子データ，初期条件，境界条件のデータ及び古い時間ステップの計算結果に基づき流体解析を実施。この処理の中で、移動物体の移動も実施する。続いてステップＳ２８００で、入力され流体解析装置記憶手段１０１に保持されたデータに基づき解析を終了するか判断する。ステップＳ２８００にて解析を終了する（ＹＥＳ）と判定した場合には、ステップＳ２９９０に進み流れ・物体移動計算を終了する。ステップＳ２８００にて、流体解析を終了しない（ＮＯ）と判定した場合は、ステップＳ２９００に進み、次に解析時間ｔを新しい時間に更新し、ステップＳ２５００に進む。ここで、図１８と図１９を用いて、ステップＳ２３００の処理内容を説明する。図１８は、２次元の解析空間に置かれた移動物体４２００と計算セル４１００，仮想粒子４３００を表示しており、図１９は、仮想粒子発生条件の入力画面を示している。図１８に示すように、仮想粒子４３００は、移動物体４２００の表面で形成される閉曲面（解析空間が２次元の空間であれば閉曲線）の内側に分布させる。仮想粒子を分布させる範囲Ｌｐを決定し、それにより決まる物体内部に決まる仮想面４４００と移動物体４２００の表面の間に、仮想粒子４３００をランダムに分布させる。ここで、発生させる仮想粒子数と仮想粒子を発生させる範囲は、入力装置１００１や外部記憶媒体１００２に保存したテキストファイルを用いて入力する。図１９は、表示装置１００３に表示された条件入力の画面４５００の例であり、少なくとも発生させる粒子の総数と粒子を発生させる範囲を指定するための数値（移動物体の表面から距離）を入力する。発生する粒子数とは、移動物体４２００の表面で形成される閉曲面内に発生させる仮想粒子の総数を指定する数字であり、発生させない場合は０と指定する。仮想粒子の数は、大きいほど精度向上が期待できる。しかし、仮想粒子数を増やすと計算負荷が大きくなるため、計算精度と計算時間の兼ね合いで決める。仮想粒子の発生範囲は、仮想粒子を発生させる範囲を指定する値であり、移動物体表面からの距離Ｌｐで指定することができる。前記Ｌｐの値は、格子点の距離オーダーとなることが必要であり、例えば、格子点間距離の平均値や前記平均値に定数倍の値を前記Ｌｐの値と決めると良い。ここでは、２次元の解析空間を用いて粒子発生処理の説明をしたが、３次元の空間に対しても同様の処理が可能である。
【００３２】
次に、図４に示したステップＳ２６００の処理について図５と図６を用いて説明する。
【００３３】
図５は、図２に示した移動物体計算手段の機能ブロック図である。移動物体計算手段２０２は、格子点区分手段２０２１，計算セル区分手段２０２２，仮想粒子区分手段２０２３，物体体積・相互作用力計算手段２０２４により構成される。格子点区分手段２０２１は、移動物体データと計算格子データに基づき、格子点４１０２が内部点であるかを判定する。判定は、光線交差法等の手法によりなされる。計算セル判定手段２０２２は、物体面データと格子点区分手段２０２１の結果に基づき、各計算セル４１００を内部セル４１０８，外部セル４１０７，境界面上セル４１０６に分類する手段である。仮想粒子区分手段２０２３は、格子点区分手段２０２３の結果と、仮想粒子データに基づき、境界面上セル４１０６に属する仮想粒子を検索する手段である。判定は、仮想粒子位置と格子点の位置から行う。また、物体体積・相互作用力計算手段２０２４は、各処理の判定結果を用いて、各計算セル４１００を占める物体の体積分率及び物体が流体へ及ぼす効果を表す相互作用力を計算する手段である。
【００３４】
図６は、図５に示した移動物体計算手段における処理を表すフローチャートである。まず、ステップＳ２６１０にて、各格子点４１０２が内部点であるかを判定する。続いてステップＳ２６２０に進み、計算セル４１００を、内部セル４１０８，外部セル４１０７，境界面上セル４１０６に区分する処理を実施する。続いて、ステップＳ２６４０にて、境界面上セル４１０６に位置する体積計算用の仮想粒子を検索する処理を実施。次に、ステップＳ２６６０に進み物体体積分率，相互作用力を計算する。
【００３５】
続いて図６に示したステップＳ２６１０からステップＳ２６４０の処理について図７から図１０を用いて説明する。
【００３６】
図７は、図６に示したステップＳ２６１０の処理を表すフローチャートである。ステップＳ２６１１にて格子点区分処理を開始し、次にステップＳ２６１２に進み、Ｍ番目（１≦Ｍ≦Ｎgm）の格子点（以下、Ｍ番格子点）の区分を開始する。続いて、ステップＳ２６１３に進み、Ｍ番格子点が内部点であるかを判定する。ステップＳ２６１３の処理は、光線交差法等を用いて実現される。続いて、ステップＳ２６１３で、物体外部にある（ＮＯ）と判定した場合、ステップＳ２６１４に進み、Ｍ番格子点の物体内部フラグをＦalseに設定し、ステップＳ２６１６に進む。ステップＳ２６１３にて内部点である（ＹＥＳ）と判定した場合、ステップＳ２６１５に進みＭ番格子点の内部フラグをＴrueに設定、続いて、ステップＳ２６１６に進む。ステップＳ２６１６では、全ての格子点の区分を完了したかを判定する。ステップＳ２６１６で、完了したと判定した場合には、ステップＳ２６１８に進み、格子点区分処理を終了する。ステップＳ２６１６にて全ての格子点の区分を完了していないと判定した場合、ステップＳ２６１７に進み、Ｍの値を更新し（例えばＭをＭ＋１の値に置き換える）、ステップＳ２６１２に進む。
【００３７】
図８は、図６に示したステップＳ２６２０の処理を表すフローチャートである。まず、ステップＳ２６２１に進み、計算セル区分処理を開始する。次にステップＳ２６２２に進み、セル番号がＫ番（１≦Ｋ≦Ｎcm）のセル（以下、Ｋ番セル）の判定を開始する。続いて、ステップＳ２６２３に進み、Ｋ番セルが外部セル４１０７であるかを判定する。セルを構成する格子点の物体内部フラグを用いることができる。セルを構成する格子点の物体内部フラグが全てＦalseであれば、外部セルである。ステップＳ２６２３にて、外部セルである（ＹＥＳ）と判定した場合、ステップＳ２６２５に進み、Ｋ番セルのセル位置フラグを外部セルと設定し、続いて、ステップＳ２６２８に進む。ステップＳ２６２３にて、外部セルでは無い（ＮＯ）と判定した場合、ステップＳ２６２４に進み、Ｋ番セルが内部セルであるかを判定する。判定には、セルを構成する格子点の物体内部フラグを用いることができる。セルを構成する格子点の物体内部フラグが全てＴrueであれば、内部セルである。ステップ２６２４で内部セルである（ＹＥＳ）と判定した場合、ステップＳ２６２６に進み、Ｋ番セルのセル位置フラグを内部セルと設定し、ステップＳ２６２８に進む。ステップＳ２６２４にて、内部セルでは無い（ＮＯ）と判定した場合、ステップＳ２６２７に進み、セル位置フラグを境界面上セルと設定、続いてステップＳ２６２８に進む。ステップＳ２６２８では、全ての計算セルを分類したかを判定し、全ての格子の分類が完了した（ＹＥＳ）と判定した場合には、ステップＳ２６３０に進み計算セル区分処理を終了する。ステップＳ２６２８にて、全ての分類が完了していない（ＮＯ）と判定した場合には、ステップＳ２６２９に進みＫの値を更新し（例えばＫをＫ＋１の値に置き換える）、ステップＳ２６２２に進む。
【００３８】
図９は、図６に示したステップＳ２６４０の処理を表すフローチャートである。最初に、ステップＳ２６４１に進み、仮想粒子区分を開始する。次にステップＳ２６４２に進み、仮想粒子番号がＪ番（１≦Ｊ≦Ｎpm）の仮想粒子（以下、Ｊ番仮想粒子）の処理を開始する。次にステップＳ２６４３に進み、Ｊ番仮想粒子の座標と計算格子データに基づきＪ番仮想粒子が所属するセル番号を検索し、その結果に基づき所属計算セル番号を設定する。続いてステップＳ２６４４に進み、Ｊ番仮想粒子が境界面上セルに位置するかを判定する。判定には、所属計算セル番号として保持している数字をセル番号として持つ計算セルのセル位置フラグを用いることができる。ステップＳ２６４４で、境界面上セルに位置しない（ＮＯ）と判定した場合、ステップＳ２６４６に進み、仮想粒子データ４３０１の境界面上フラグをＦalseと設定し、ステップＳ２６４８に進む。ステップＳ２６４４にて境界面上セルに位置する（ＹＥＳ）と判定した場合、ステップＳ２６４５に進み、境界面上セル内の仮想粒子数を保存する変数であるＮpcを更新し、ステップＳ２６４７に進み、境界面上フラグをＴrueに設定する。続いて、ステップＳ２６４８に進む。ステップＳ２６４８では、全ての仮想粒子の分類を完了したかを判定し、全ての仮想粒子の分類を完了した（ＹＥＳ）と判定した場合には、ステップＳ２６５０に進み、仮想粒子検索処理を終了する。ステップＳ２６４８にて、全ての仮想粒子の分類を完了していない（ＮＯ）と判定した場合には、ステップＳ２６４９に進みＪの値を更新し（例えばＪの値をＪ＋１に置き換える、ステップＳ２６４２に進む。
【００３９】
図１０は、図６に示したステップＳ２６６０の処理を示すフローチャートである。最初に、ステップＳ２６６１に進み、体積分率・相互作用力の計算を開始する。次に、ステップＳ２６６２に進み、セル番号がＩ番（１≦Ｉ≦Ｎcm）である計算セル（以下、Ｉ番セル）の計算を開始する。次にステップＳ２６６３に進み、Ｉ番セルが内部セルかを判定する。この判定には、セルデータ４１０５のセル位置フラグを用いることができる。ステップＳ２６６３にて、内部セルである（ＹＥＳ）と判定した場合、ステップＳ２６６５に進み、物体体積分率を１に設定、前記物体移動速度ｖeｃ(ｕP)と前記式（３）を用いて相互作用力を計算し、ステップＳ２６６８に進む。ステップＳ２６６３にて、内部セルでは無い（ＮＯ）と判定した場合、ステップＳ２６６４に進み、Ｉ番セルが外部セルであるかを判定する。ステップＳ２６６４にて、物体外部セルである（ＹＥＳ）と判定した場合には、ステップＳ２６６７に進み、物体体積分率と相互作用力を０と設定し、ステップＳ２６６８に進む。ステップＳ２６６４にて、物体外部セルでは無い（ＮＯ）と判定した場合、ステップＳ２６６６進み、仮想粒子データに基づく体積分率と相互作用力の計算を行う。以下に、ステップＳ２６６６における処理の具体例を述べる。境界面上セルに区分される計算セルに所属する仮想粒子の総数をＮST、Ｉ番セルに含まれる仮想粒子個数をＮPI、移動物体の総体積をＶm、Ｉ番セルの体積をＶI、内部セルに区分された計算セルの体積の総和をＶIN、ｊ番仮想粒子の速度をｖeｃ(ｕPj)、移動物体の移動速度をｖeｃ(ｕP)、流体の密度をρとする。本定義のもと、Ｉ番セルの体積分率θI，相互作用力ｖeｃ(ＦSI)を式（５）から式（７）で求める。
【００４０】
【数５】
2011070333.tif 000007
【００４１】
【数６】
2011070333.tif 000008
【００４２】
【数７】
2011070333.tif 000009
【００４３】
ここで、Δｔはシミュレーションにおける時間刻み幅であり、計算の安定性や求める計算精度に応じて決定する。式（６）の総和は、Ｉ番セルに所属する粒子について計算する。
【００４４】
式（５）及び式（６）では、仮想粒子の影響が及ぶ範囲を点に限定しているが、分布関数を用いることで、仮想粒子の影響が及ぶ範囲を広げることができ、さらに計算精度の向上が期待できる。例えば、式（８）から式（１０）のように計算ができる。
【００４５】
【数８】
2011070333.tif 000010
【００４６】
【数９】
2011070333.tif 000011
【００４７】
【数１０】
2011070333.tif 000012
【００４８】
ここで、ｖeｃ(ｘ)及びｖeｃ(ｘp)は解析空間中の位置を表すベクトル、ｖeｃ(ｘp,j)は解析空間中の仮想粒子の位置を表すベクトル、σは粒子を分布させる範囲の定数倍で与える実数、である。また、式（８）と式（９）の総和は、Ｉ番セルに所属する粒子について計算する。以上に示したのは、一例であり、ステップＳ２６６６の体積分率の計算には、これらの方法に限ることなく任意の方法を選択できる。
【００４９】
図２０は、本発明の実施形態である流体解析装置を用いて計算した体積分率の表示画像４６００であり、２次元解析空間４７００における計算セル４１００，移動物体４２００，仮想粒子４３００及び各計算セルにおける体積分率計算結果を示している。本発明により、外部セル４１０７で体積分率を０、内部セル４１０８で体積分率を１、境界面上セル４１０６で体積分率を０より大きくかつ１より小さく計算することができるため、外部セル４１０７における物体体積分率が０、内部セル４１０８における物体体積分率が１、境界面上セル４１０６における物体体積分率が０より大きく１より小さい、と表示できる。また、移動物体４２００の表面内部かつ移動物体４２００の表面近傍に仮想粒子４３００が分布する結果を画面上に表示することができる。
【符号の説明】
【００５０】
１００  入力データ記憶手段
１０１  流体解析条件記憶手段
１０２  計算格子データ記憶手段
１０３  移動物体データ記憶手段
１０４  仮想粒子データ記憶手段
２００  流れ・物体移動解析手段
２０１  流れ場解析手段
２０２  移動物体計算手段
２０３  仮想粒子発生手段
３００  結果表示手段
１０００  流体解析装置
１００１  入力装置
１００２  外部記憶媒体
１００３  表示装置
１００４  解析装置
１００５  主記憶装置
１００６  中央処理装置
２０２１  格子点区分手段
２０２２  計算セル区分手段
２０２３  仮想粒子区分手段
２０２４  体積分率・相互作用力計算手段
４１００  計算セル
４１０１  セル表面
４１０３  格子点データ
４１０４  セル表面データ
４１０５  セルデータ
４１０６  境界面上セル
４１０７  外部セル
４１０８  内部セル
４１１１  内部点
４１１２  外部点
４２００  移動物体の表面
４２０１  物体微小面
４２０２  物体点
４２０３  物体点データ
４２０４  物体微小面データ
４３００  仮想粒子
４３０１  仮想粒子データ
４４００  仮想粒子分布範囲
４５００  入力画面画像
４６００  出力画面画像
【図１】
2011070333.tif 000013
【図２】
2011070333.tif 000014
【図３】
2011070333.tif 000015
【図４】
2011070333.tif 000016
【図５】
2011070333.tif 000017
【図６】
2011070333.tif 000018
【図７】
2011070333.tif 000019
【図８】
2011070333.tif 000020
【図９】
2011070333.tif 000021
【図１０】
2011070333.tif 000022
【図１１】
2011070333.tif 000023
【図１２】
2011070333.tif 000024
【図１３】
2011070333.tif 000025
【図１４】
2011070333.tif 000026
【図１５】
2011070333.tif 000027
【図１６】
2011070333.tif 000028
【図１７】
2011070333.tif 000029
【図１８】
2011070333.tif 000030
【図１９】
2011070333.tif 000031
【図２０】
2011070333.tif 000032