JP5033211B2 - 流体シミュレーションにおける境界位置決定方法 - Google Patents

Description

本発明は、流体シミュレーションにおける境界位置決定方法に関し、詳しくは流体モデルと物体モデルとの境界位置を高速に設定することができる技術に関する。

近年、様々な場面で、コンピュータを用いた流体シミュレーションが行われている。流体シミュレーションの一例として、ゴルフボールの周囲の空気の流れを解析するものや、ゴム押出機の内部のゴムの流れを解析するもの等が知られている。このような流体シミュレーションは、飛行特性に優れたゴルフボールのディンプルの開発や、押出し抵抗の少ない押出機ないしゴム配合の開発等に役に立つ。

流体シミュレーションにおいては、通常、三次元の固体を有限個の要素で表した物体モデルと、該物体モデルの周囲を流れる流体の空間としてメッシュの流体解析用の格子点からなる流体モデルとが定義される。そして、前記各格子点において、流体の温度、圧力及び速度などが計算される。関連する文献としては次のものがある。

特開２００４−２９９４３３号公報

ところで、流体シミュレーションでは、通常、物体モデルの節点と流体モデルの格子点とは異なる配列を持っている。従って、流体モデルのどの格子点が物体モデルの内部にあり、どの格子点が物体モデルの外部にあるのかを知る必要がある。つまり、流体シミュレーションを行うに際しては、それに先立ち、流体モデル側に、物体モデルとの境界点の位置を定義するプロセスが必要になる。

従来、このような境界を流体モデル側に定義するために、物体モデルの表面の各要素の節点に最も近い格子点を探索するいわゆる”点対点”のサーチが行われており、代表的なサーチ手法として、最近傍探索法が知られている。最近傍探索の高速化手法は、最近傍候補の絞り込みによるアプローチと、距離計算の打ち切りによるアプローチの２種類に大別される。

しかしながら、最近傍探索法では、流体モデルの大型化により探索する対象が増えると、最近傍点の探索に時間がかかるという問題があった。また、近年では、ＡＮＮ（Approximate Nearest Neighbor）やｋｄ−Ｔｒｅｅといったアルゴリズムも提案されているが、やはり格子点の数が増えると、最近傍点の探索に多くの時間が必要である。

さらに、物体モデルが移動及び／又は変形するような非定常計算では、各ステップで、流体モデルと物体モデルとの境界点を逐次再認識する必要がある。従って、効率の悪い境界認識処理では、計算コストが莫大になるという問題があった。

本発明は、以上のような問題点に鑑み、案出なされたもので、流体モデルに、物体モデルとの境界位置を能率良くかつ高速に定義することができる流体シミュレーションの境界位置決定方法を提供することを主たる目的としている。

本発明のうち請求項１記載の発明は、三次元の固体を有限個の要素で表した物体モデルと、流体が流れる空間を規定する流体モデルとを用いて少なくとも前記流体の流れを解析する流体シミュレーションを行うに際し、
ａ）コンピュータに、表面が複数の平面からなる前記物体モデルを入力するステップ、
ｂ）前記コンピュータに、前記物体モデルを覆うように三次元空間に多数の流体解析用の格子点が定義された前記流体モデルを入力するステップ、
ｃ）前記コンピュータに、前記流体モデルの前記格子点を通りかつ前記流体モデルを貫く複数の直線を定義するステップ、
ｄ）前記コンピュータが、前記物体モデルの表面と前記直線との交点の座標値を計算しかつこの交点の座標値を記憶するステップ、及び
ｅ）前記コンピュータが、上記ステップｄで得られた交点の座標値に基づいて、前記流体モデルの全ての格子点を固体領域又は流体領域のいずれかに決定するステップ
を含んで前記流体モデルと前記物体モデルとの境界位置を前記流体モデルに定義し、
前記ステップｅは、前記コンピュータの記憶部に、探索用データベースを作る前処理ステップと、
前記コンピュータが前記探索用データベースを使用して前記交点の座標値に最も近い前記格子点の座標値の探索を行う探索ステップとを含み、
前記探索用データベースは、前記直線上の格子点の座標値がなす実数群の区間（x[1]〜x[n]）を小領域であるバケットで隙間無くかつ等間隔でＱ個に区分されたバケット列を有し、
前記各バケットには、先頭バケットから末尾バケットまで順番に割り当てられた固有のバケット番号と、
各バケットの中に含まれる前記実数の個数であるバケットサイズとが記憶されるとともに、
少なくとも前記先頭バケットを除いた各バケットには、前記バケットサイズがゼロでなくかつ当該バケットに最も近いバケットのバケット番号を示すラストフィルドバケットが記憶され、
前記探索ステップは、前記交点の座標値が属するバケットを計算する第１のステップと、
前記計算されたバケットのバケットサイズを参照する第２のステップと、
前記参照されたバケットサイズが０でない場合、当該参照されたバケットから前記交点の座標値の最近傍値を計算する一方、
前記参照されたバケットサイズが０の場合、そのバケットのラストフィルドバケットが示すバケットから前記交点の座標値の最近傍値を計算する第３のステップとを含むことを特徴とする。

また請求項２記載の発明は、前記流体モデルの前記格子点は、ｘ−ｙ−ｚの三次元空間のｘ軸方向をｓ分割、ｙ軸方向をｍ分割、ｚ軸方向をｎ分割したメッシュの交差点上に設けられ、前記直線は、前記メッシュに一致して定義される請求項１に記載の流体シミュレーションにおける境界位置決定方法である。

また請求項３記載の発明は、前記流体モデルの格子点は、ｒ−θ−ｚの円筒形空間のｒ方向をｓ分割、θ方向をｍ分割、ｚ方向をｎ分割したメッシュの交差点上に設けられ、前記直線は、前記メッシュに一致して定義される請求項１記載の流体シミュレーションにおける境界位置決定方法である。

また請求項４記載の発明は、前記流体モデルの格子点は、ｒ−θ−φの球形空間のｒ方向をｓ分割、θ方向をｍ分割、φ方向をｎ分割したメッシュの交差点上に設けられ、前記直線は、前記メッシュに一致して定義される請求項１記載の流体シミュレーションにおける境界位置決定方法である。

また請求項５記載の発明は、前記最近傍値は、前記交点の座標値に最も近い物体モデル内の格子点の座標値とし、前記探索用データベースには、前記ラストフィルドバケットとして、当該バケットよりも先頭バケット側で当該バケット番号に最も近いバケット番号が記憶されている請求項１記載の流体シミュレーションにおける境界位置決定方法である。

また請求項６記載の発明は、前記探索用データベースの少なくともバケットサイズが０でないバケットには、そのバケット内の実数のうちの最小値及び最大値がさらに記憶されている請求項１乃至５のいずれかに記載の流体シミュレーションにおける境界位置決定方法である。

また請求項７記載の発明は、前記第３のステップにおいて、参照されたバケットサイズが０でない場合、前記交点の座標値と、そのバケットの最小値とを比較し、前記交点の座標値が最小値よりも小さい場合には、そのバケットのラストフィルドバケットの最大値を最近傍値とする請求項６記載の流体シミュレーションにおける境界位置決定方法である。

また請求項８記載の発明は、前記第３のステップにおいて、参照されたバケットサイズが０でない場合、前記交点の座標値と、そのバケットの最小値とを比較し、前記交点の座標値が最小値以上の場合には、そのバケットに含まれる実数に対してx[ ipoint ] ≦ｑ＜x[ ipoint + 1 ]となるx[ ipoint ] を最近傍値とする請求項６又は７記載の流体シミュレーションにおける境界位置決定方法である。

本発明によれば、流体モデルに定義された格子点を通る直線と、物体モデルの表面をなす平面の要素との交点に基づいて、流体モデルの全ての格子点が固体領域又は流体領域のいずれかに決定される一次元のサーチアルゴリズムが提供される。従って、従来のように、物体モデルの表面の要素の節点に最も近い格子点を探す点対点の多次元のサーチアルゴリズムに比べて次元数が減り、計算が大幅に簡素化できる。従って、流体モデルと物体モデルとの境界を認識するための計算時間が大幅に短縮され、高速な処理が可能になる。また、境界となる前記直線上の交点間に含まれる格子点については、一括して非流体領域（流体モデルと物体モデルの交差部分）として決定できるため、格子点に非流体領域を定義する処理であるいわゆるソリッドマーキング処理についても大幅に高速化できる。

本実施形態のシミュレーション方法を実施するためのコンピュータ装置の構成図である。 本実施形態の流体シミュレーション方法の処理手順の一例を示すフローチャートである。 物体モデルの一例を視覚化して示す部分斜視図である。 流体モデルの一例を視覚化して示す斜視図である。 （ａ）は図４の側面図、（ｂ）はそのＡ部拡大図である。 境界認識処理の概要を示す斜視図である。 境界認識処理の処理手順の一例を示すフローチャートである。 物体モデルと要素平面との交点を説明する線図である。 （ａ）、（ｂ）は物体モデルと直線の交差を説明する線図である。 物体モデルと流体モデルとの境界を説明する線図である。 物体モデルと直線との交差を説明する線図である。 探索方法の処理手順を示すフローチャートである。 実数群及び探索用データベースの構造を示す線図である。 探索用データベースを作成する処理手順を示すフローチャートである。 ラストフィルドバケットを決定する処理手順の一例を示すフローチャートである。 探索ステップの処理手順の一例を示すフローチャートである。 円筒形座標の流体モデルの実施形態を示す線図である。 球形座標の流体モデルの実施形態を示す線図である。

以下、本発明の実施の一形態が図面に基づき説明される。
図１には、本実施形態の流体シミュレーション方法を実施するためのコンピュータ装置１が示されている。このコンピュータ装置１は、本体１ａと、入力手段としてのキーボード１ｂ、マウス１ｃと、出力手段としてのディスプレイ装置１ｄとを含んで構成されている。本体１ａには、演算処理装置（ＣＰＵ）、作業用メモリ及び磁気ディスク等の補助記憶手段の他、ＣＤ−ＲＯＭやフレキシブルディスクのドライブ１ａ１、１ａ２などを適宜具えている。そして、前記磁気ディスクには本発明に係る処理を実行するためのプログラムが記憶されている。

図２には、本発明の流体シミュレーションにおける境界位置決定方法の処理手順の一例が示される。以下、順に説明する。

本実施形態の流体シミュレーション方法では、先ず、上記コンピュータ装置１に、物体モデルが入力される（ステップＳ１）。

前記物体モデルは、固体（非流体物）を有限個の要素でモデル化したものである。図３には、本実施形態の物体モデル２を視覚化した斜視図が示される。該物体モデル２は、三次元の物体としてのゴルフボールをモデル化したゴルフボールモデル２ａからなる。

本実施形態の流体シミュレーションでは、ゴルフボールモデル２ａの周囲の流体（空気）の流れのみが観察される。従って、本実施形態のゴルフボールモデル２ａは、多数の平面要素ｅを連ねて球形に形成された殻構造で十分であり、内部構造は省略されている。これは、モデル２ａの作成時間を短縮するのに役立つ。また、変形が考慮されないため、ゴルフボールモデル２ａは剛表面で作られても良い。ただし、何らかの変形計算も合わせて行うような場合には、ゴルフボールモデル２ａは、例えば有限要素法等の数値解析法で計算が可能なように、例えば変位を未知数とするラグランジェ（Lagrange）要素が好適に用いられる。

また、コンピュータ装置１に物体モデル２を入力するとは、少なくともゴルフボールモデル２ａの表面の形状データがコンピュータ装置１の記憶手段に数値データとして入力されかつ記憶されることを意味する。この作業は、前記入力手段を用いてオペレータにより行われる。

次に、流体モデル３がコンピュータ装置１に入力される（ステップＳ２）。
流体モデル３は、図４、図５に示されるように、流体としての空気が流れる三次元の空間領域を定める。本実施形態の流体モデル３は、オイラー（Euler）要素で構成される。オイラー要素とは、図５（ａ）のＡ部拡大図である図５（ｂ）に示されるように、前記三次元の空間に予め定められた流体解析用の格子点３ａを有し、各格子点３ａでの空気の圧力、温度及び速度などが未知数とされる。

本実施形態の流体モデル３は、ゴルフボールモデル２ａを完全に覆う大きさを有した略直方体状に定義されている。換言すれば、流体モデル３は、物体モデル２とオーバーラップするように設けられる、これにより、必要な格子点３ａの物理量を計算することで、ゴルフボールモデル２ａの周囲に流れる空気の状態を知るのに役立つ。ただし、流体モデル３は、観測したい部位に合わせて種々その大きさを設定できるのは言うまでもない。

本実施形態の流体モデル３の格子点３ａは、ｘ−ｙ−ｚの三次元空間のｘ軸方向をｓ分割、ｙ軸方向をｍ分割、ｚ軸方向をｎ分割した直交格子のメッシュＭの交差点上に設けられている。ｓ、ｍ及びｎは、いずれも自然数であり、これらはモデルの大きさ、解析の精度、流速等に応じて適宜定められ、あくまで一例であるが、ゴルフボールモデル２ａの外径が４２．７mmの例では、各軸において、格子点のピッチは、例えば２．５×１０^-2〜１mm程度が望ましい。

次に、本実施形態では、流体モデル３に、該流体モデル３と前記ゴルフボールモデル２ａとの境界を定義する境界認識処理が行われる（ステップＳ３）。

図５（ｂ）のハッチングにて示されるように、流体モデル３のゴルフボールモデル２ａとオーバーラップする格子点３ａには、空気は流れない。従って、上で述べたように、流体シミュレーションを行うに際して、流体モデル３のどの格子点３ａが物体モデル２の内部（流体が流れない固体領域）に位置し、どの格子点３ａが物体モデル２の外部（流体が流れる流体領域）にあるのかを決定する必要がある。そのために、流体モデル３に、境界点が決定される。

本実施形態の境界認識処理では、図６に示されるように、概略、流体モデル３に格子点３ａを通る多数の直線Ｌを引き、これらの直線Ｌとゴルフボールモデル２ａの表面との交点Ｃを求め、これに基づいて、前記境界が定義される。以下、このような処理の具体的な実施形態が説明される。

図７には、コンピュータ装置１で行われる境界認識処理の一例としてのフローチャートが示される。本実施形態では、先ず、前記流体モデル３上に、多数の直線Ｌが定義される（ステップＳ３１）。前記直線Ｌは、流体モデル３の格子点３ａを通りかつ流体モデル３を貫く長さで定義される。ここで、「貫く」とは、直線Ｌが各ｘ、ｙ、ｚ軸方向の流体モデル３の大きさ以上であれば足りるものとする。

本実施形態では、前記直線Ｌは、前記メッシュＭに一致して定義される。具体的には、図６に示したように、ｘ軸、ｙ軸及びｚ軸と平行な３種類Ｌ１、Ｌ２及びＬ３を含み、これらは前記各軸方向の格子点の分割数に一致して定義される。即ち、本実施形態のような直交座標系では、各ｘ、ｙ及びｚ軸に沿って３方向に直線が配され、その配置及び数は、次のようになる。
１）ｓ×ｍ個の直線がｚ軸と平行に配置
２）ｓ×ｎ個の直線がｙ軸と平行に配置
３）ｍ×ｎ個の直線がｘ軸に平行に配置
従って、本実施形態では、全ての格子点３ａは、いずれかの直線Ｌ上に位置することになる。

次に、変数ｉが１にセットされ（ステップＳ３２）、物体モデル２（ゴルフボールモデル２ａ）の表面のｉ番目の要素平面の節点の座標が作業用メモリに読み込まれる（ステップＳ３３）。

図８には、ゴルフボールモデル２ａの表面に位置する三角形平面の要素ｅ…の一つが抜き出して描かれている。このステップＳ３３では、前記要素ｅの平面（要素平面ｎ）の３つの節点Ｐ１、Ｐ２及びＰ３の各座標が作業メモリに読み込まれる。なお、コンピュータ装置１は、この処理を迅速に行うために、要素平面ｎの節点座標が記憶されたデータファイルを予め具えておくことが望ましい。

次に、コンピュータ装置１は、全ての直線Ｌを特定するための情報（方程式等）が、例えば作業用メモリに読み込まれる（ステップＳ３４）。

次に、各直線Ｌとゴルフボールモデル２ａの前記要素平面ｎとが平行か否かが判断される（ステップＳ３５）。直線Ｌと要素平面ｎとが平行の場合（ステップＳ３５でＹ）、直線と要素平面ｎとは交差しないものとして判断される（ステップＳ３７）。

他方、ステップＳ３５でＮ、即ち、直線Ｌが要素平面ｎと平行でない場合、直線Ｌと要素平面ｎを含む平面Ｎとの交点の座標値が計算される（ステップＳ３６）。図８に示したように、要素平面ｎを含む平面Ｎは、前記３つの節点の座標Ｐ１乃至Ｐ３を用いた方程式により特定される。また、この平面Ｎと、直線Ｌとの交点も、簡単な計算で求めることができるのは言うまでもない。

次に、ステップＳ３６で得られた交点が、要素平面ｎ内にあるか否かが判断される（ステップＳ３８）。そして、交点Ｃが要素平面ｎ内にある場合（ステップＳ３８でＹ）、その交点Ｃはゴルフボールモデル２ａの表面の交点として記憶される（ステップＳ３９）。

そして、変数ｉを１ずつ増加させながら、上記の処理が、全ての直線Ｌと全ての要素平面ｎとの組合せについて行われる（ステップＳ４０）。

以上の処理により、ゴルフボールモデル２ａの表面を構成している要素平面ｎの全てについて、ステップＳ３１で定義された直線が交差するかどうかが探索される。また、全ての直線Ｌについて、ゴルフボールモデル２ａの表面を構成している要素平面と交差したかどうかがチェックされる。なお、ステップＳ３９の処理において、直線Ｌが複数回ゴルフボールモデル２ａの表面と交差した場合には、どの順番で交差したかについても記憶される。

次に、上記ステップＳ３９で得られた物体モデル２の表面の交点Ｃの座標値に基づいて、流体モデル３に、該流体モデル３と前記物体モデル２との境界が定義される（ステップＳ４１）。

本実施形態の境界の定義は、図９（ａ）、（ｂ）に示されるように、三次元空間をのびる無限長と定義できる直線Ｌは、物体モデル２（固体部分）を横切るか、又は全く横切らないか、のいずれかであるという事実に基づいている。また、図９（ａ）に示されるように、直線Ｌが物体モデル２を横切るとき、直線Ｌは、物体モデル２の表面と少なくともＡ、Ｂの２点で交差する。

従って、図１０に代表的な例が示されるように、コンピュータ装置１は、１本の直線Ｌｉ上に交点Ｃが２つのみである場合、前記直線Ｌｉにおいて、２つの交点Ｃ、Ｃ間に含まれる全ての格子点３ａを流体が流れない固体領域Ｓ（黒塗り丸）として定義しうる一方、上記以外の格子点３ａを流体が流れる流体領域Ｆ（白抜き丸）として定義することができる。そして、流体領域Ｆのうち、物体モデル２に最も近い格子点３ａが計算上の境界点となる。なお、物体モデル２の表面上に格子点３ａが一致する場合、その格子点３ａは固体領域として取り扱われる。

また、図１１に示されるように、物体モデル２の形状によっては、直線Ｌは、２回以上、物体モデル２の表面と交差する場合もある。この場合、交差する回数は必ず偶数回となり、図１１の例では１４回である。従って、任意の方向から数えて奇数番目の交点は、物体モデル２の表面から物体モデル２の内部へ貫通する開始点とみなすことができる。逆に、偶数番目の交点は、物体モデル２の貫通を終える終了点とみなすことができる。

従って、コンピュータ装置１は、前記交点がｋ個以上（ただし、ｋは４以上の偶数）ある直線の場合、最初に物体モデル２に交差した交点から順番に各交点をＣ1、…、Ｃkとするとき、交点ＣiからＣi+1（ただし、ｉは１から始まるｋ−１までの奇数である）の間に含まれる格子点３ａを流体が流れない固体領域Ｓとし、かつ それ以外の格子点３ａを流体が流れる流体領域Ｆとして定義することができる。

以上のように、直線Ｌ上の交点について、その座標と、物体モデル２への交差順番とが分かれば、形状に拘わらず、物体モデル２の内部である固体領域Ｓを簡単に特定することができる。このように、本発明の境界認識処理によれば、流体モデル３に定義された直線Ｌと、物体モデル２の表面をなす要素平面ｎとの交点Ｃを用い、この交点Ｃの座標値を用いて境界が定義される。従って、従来のように、物体モデル２の表面の要素の節点に最も近い格子点を探す多次元上での点対点のサーチアルゴリズムに比べて次元数が減り、計算が大幅に簡素化できる。従って、計算時間が大幅に短縮される。また、本実施形態によれば、境界となる前記直線Ｌ上の交点間に含まれる格子点３ａについては、一括して非流体領域Ｓとして決定できる。従って、格子点に非流体領域Ｓを定義していくソリッドマーキング処理についても大幅に高速化しうる。

また、物体モデル２が移動又は回転した場合には、上記境界認識処理が繰り返され、新たな境界が流体モデル３に定義される。

次に、前記ステップＳ４１の処理を高速に行う好ましい実施例について述べる。
この実施形態では、実数が小さいものから順番に並んでいる実数群（x[1],x[2], x[3] … x[n]）に対して、ある探索対象となる探索点ｑが与えられ、その探索点ｑに最も近い最近傍値を前記実数群の中から探索する方法が利用される。即ち、上記実数群に各直線Ｌ上に並ぶ格子点３ａ…の一次元の座標値が用いられる一方、探索点ｑには前記交点Ｃの座標値が用いられる。

この実施形態の探索方法では、図１２に示されるように、探索用データベースを作る前処理ステップＳ４１ａと、前記コンピュータ装置１が前記探索用データベースを使用して前記探索を行う探索ステップＳ４１ｂとが含まれる。

前記前処理ステップＳ４１ａは、探索ステップＳ４１ｂの前に一度だけ実行されれば良い。そして、この前処理ステップＳ４１ａにおいて、探索ステップＳ４１ｂで必要になる探索用データベースがコンピュータ装置１に作成される。

探索用データベースは、図１３に視覚化して示される。図１３では、前記直線Ｌ上の格子点がなす一次元の座標が小さいものから順番に並んでいる実数群（x[1],x[2], x[3] … x[n]）がある。前記探索用データベースは、前記実数群の区間（x[1]〜x[n]）を小領域であるバケットＢで隙間無くかつ等間隔でＱ個に区分されたバケット列Ｂｒを有する。

［バケット］
バケットＢは、１次元の空間と言え、２つの実数ｙ1及び2を用いてｙ1≦ x[] ＜ｙ2を満たす空間として表すことができる。各バケットはいずれも同じ空間サイズを有し、前記実数群の区間（x[1]〜x[n]）を隙間無くかつ等間隔で分割している。

［バケット番号］
また、各バケットＢには、先頭バケットから末尾バケットまで順番に割り当てられた固有のバケット番号が定義される。この実施形態では、先頭バケットのバケット番号が１、末尾バケットのバケット番号がＱ、これらの間でバケット番号が１ずつ増加している。

［バケットサイズ］
さらに、各バケットには、その中に含まれる前記実数の個数であるバケットサイズが記憶される。この実施形態では、バケット番号１、２、６、９及び１０のバケットにそれぞれ１つの実数x[1]、x[2]、x[3] 、x[4] 、x[5] が含まれているので、これらのバケットサイズにはいずれも”１”が記憶される。また、実数を一つも含まないバケットには、”０”が記憶される。

［ラストフィルドバケット］
さらに、少なくとも前記先頭バケットを除いた各バケットには、ラストフィルドバケットが記憶される。ラストフィルドバケットは、前記バケットサイズがゼロでなくかつ当該バケット番号に最も近い近接バケット番号を示すラストフィルドバケットが記憶される。ここで「最も近い」とは、３つの意味を含むものとする。一番目の意味は、文字通り、当該バケットに最も近いという意味、二番目の意味は、当該バケットの先頭バケット側で最も近いという意味、三番目の意味は当該バケットの末尾バケット側で最も近いという意味である。この方法では、これらのいずれかの意味を選択して、実施することができる。

このラストフィルドバケットは、コンピュータ装置１が、実数群をx[1]からx[n]までスキャンし、ｙ1≦ x[i] ＜ｙ2を満たすか否かを調べていく課程の中で、最後にバケットサイズが≠０となったときのバケットを、メモリ上に「ラストフィルドバケット」として記憶させておく。なお先頭バケットには、ラストフィルドバケットとして、自らのバケット番号を記憶させることができる。

さらに、探索用データベースの少なくともバケットサイズがゼロでないバケットＢには、そのバケット内の実数のうちの最小値及び最大値がさらに記憶されている。バケットサイズが１の場合、最小値及び最大値はいずれも同じ値である。

次に、この前処理ステップＳ１において探索用データベースを作成する具体的な処理手順の一例が図１４に示される。

先ず、コンピュータ装置１は、前記実数群（x[1],x[2], x[3] … x[n]）をメモリ上に読み込む（ステップＳ５１）。

次に、コンピュータ装置１は、前記実数群の区間（x[1]〜x[n]）を前記バケットで隙間無くかつ等間隔でＱ個に区分する（ステップＳ５２）。これにより、１、２、３、… Ｑのバケット番号を有するバケットが作られる。前記Ｑは、０よりも大きい整数であり、ユーザにより決定される。そして、整数Ｑは、予め入力手段５からコンピュータ装置１に入力されている。通常、前記整数Ｑは、実数の組（x[1]、x[2]、 x[3] … x[n])の個数ｎの数倍大きな値をとるように定められるのが望ましい。

次に、コンピュータ装置１は、入力された前記整数Ｑに基づき、前記メモリ上に、バケット番号、バケットサイズ、ラストフィルドバケット、最小値及び最大値を記憶するための領域をそれぞれＱ個割り当てる（ステップＳ５３）。

次に、コンピュータ装置１は、前記実数群（x[1],x[2], x[3] … x[n]）を各々のバケットに代入するとともに、割り当てられたメモリの各領域にそれぞれの値を記憶させる（ステップＳ５４）。この処理では、先ず、バケット総数を示す前記整数Ｑと、実数の最小値x[1]、同最大値x[n]とを用いて、下式（１）にて各バケットの間隔Δxが計算される。
Δx＝（x[n] − x[1]）／Ｑ …（１）

次に、コンピュータ装置１は、実数x[1]からx[n]まで、順番に次の計算を行う。即ち、現在の実数の値をx[i]とすると、x[i]がどのバケットに属しているのかを見つけるため、下記式（２）の計算を行う。
ibucket ＝ ＃切捨て(（x[i] − x[1]）／ Δx)＋１ …（２）
（ここで、「＃切り捨て」は、関数であり、戻り値として、小数点以下を切り捨てた整数値を返すものとする。）

ここで、 ibucket のバケット番号を持つバケット（これを[ibucket]として示す。）には、１つの実数x[i] を含むことが分かる。従って、コンピュータ装置１は、下記式（３）により、当該バケット[ibucket]に記憶されているバケットサイズの値を一つ増やす処理を行う。
バケットサイズ[ibucket] = 現在のバケットサイズ＋１ … （３）

また、この処理で新たにバケット[ibucket]内にx[i]が入るので、中央演算装置２は、最小値又は最大値を書き換える必要が生じたときには、それらの値を書き換える。

次に、バケットサイズ、最小値、最大値が全てのバケットに対して決定されるが、この例では、いくつかのバケットは、そのバケットサイズが０となり、実数x[]とは無関係になっていることが分かる。この様なバケットサイズ＝０のバケットについては、無視して検索対象から除外するのが得策である。

次に、コンピュータ装置１は、先頭バケットを除いて、各バケットのラストフィルドバケットを計算し、かつ、記憶させる（ステップＳ５５）。この実施形態では、ラストフィルドバケットは、バケットサイズがゼロでなくかつ当該バケットに先頭バケット側で最も近い（即ち、前記二番目の意味）バケットのバケット番号として定義される。これは、現在のバケットのバケット番号をｋ、ラストフィルドバケットをｐとすると次の条件ａないしｃで表すことができる。
（ａ）ｐ＜ｋ
（ｂ）バケット番号ｐのバケットサイズ≠０
（ｃ）ｋ−ｐが最小値

また、コンピュータ装置１がラストフィルドバケットを決定する処理手順の一例は図１５に示される。

コンピュータ装置１は、変数ＬＦＢ及び変数ｋにそれぞれ１を代入し（ステップＳ６０、Ｓ６１）、バケット番号ｋについて、バケットサイズが０より大きいかを判断する（ステップＳ６２）。

ステップＳ６２で結果が真（Ｙ）の場合、コンピュータ装置１は、バケット番号ｋのラストフィルドバケットの値に、現在の変数ＬＦＢの値を代入する（ステップＳ６３）。その後、コンピュータ装置１は、変数ＬＦＢの値を、現在のバケット番号ｋに更新し（ステップＳ６４）、ステップＳ６５を実行する。

他方、ステップＳ６２で結果が偽（Ｎｏ）の場合、バケット番号ｋのラストフィルドバケットの値に、現在の変数ＬＦＢの値を代入するが（ステップＳ６６）、変数ＬＦＢの値は更新せずに、ステップＳ６５を実行する。

ステップＳ６５では、現在のバケット番号ｋが末尾バケット（即ち、ｋ＝Ｑ）かを判断し、結果が真（Ｙｅｓ）の場合には処理を終える。他方、ステップＳ６５の結果が偽（Ｎｏ）の場合、変数ｋに１を加算し（ステップＳ６７）、末尾バケットになるまでステップＳ２２以降を繰り返す。

以上の処理により、前処理ステップが完了し、探索用データベースが作成される。なお、探索用データベースは、アクセス速度が早いメモリ上に設けられるのが好ましいが、ハードディスク等の補助記憶装置に作成されても良い。

前記探索ステップでは、前記探索用データベースを使用して、図１６の処理が行われる。先ず、１本の直線Ｌに関して、探索点ｑとしての交点の座標値の一つが与えられたとき（ステップＳ７０）、コンピュータ装置１は、前記探索点ｑが属するバケットのバケット番号（ibucket）を計算する（第１のステップ：ステップＳ７１）。この計算は、バケットが等間隔で設定されているため、単純な下記式（４）の一次式で所定のバケット見つけることができる。
ibucket ＝ ＃切捨て(（ｑ − x[1]）／ Δx)＋１ …（４）
ここで、上記「＃切捨て」は、関数であり、戻り値として、小数点以下を切り捨てた整数値を返すものとする。また、Δｘは、上記式（１）で表される。

次に、コンピュータ装置１は、前記計算されたバケット（ibucket）のバケットサイズを参照し、その値が０か否かを判断する（ステップＳ７２）。もし、バケットサイズが０でない場合（ステップＳ７２でＮ）、当該参照されたバケットの中から最近傍値を計算する（ステップＳ７３、Ｓ７４）。

具体的には、参照されたバケットサイズが０でない場合、前記探索点ｑと、そのバケットの最小値とを比較し、探索点ｑが最小値よりも小さい場合（ステップＳ７３でＹｅｓ）、そのバケットのラストフィルドバケットの最大値を最近傍値とする（ステップＳ７５）。

また、ステップＳ７３でＮ、即ち、バケットサイズが０でなく、かつ、探索点ｑがそのバケットの最小値以上の場合には、そのバケットに含まれる実数に対してx[ ipoint ] ≦ｑ＜x[ ipoint + 1 ]となるx[ ipoint ] を最近傍値とする（ステップＳ７４）。

一方、中央演算装置２は、前記参照されたバケットサイズが０の場合（ステップＳ７２でＹｅｓ）、そのバケットのラストフィルドバケットが示すバケットから前記最近傍値を計算する（ステップＳ７５）。

このように、本実施形態の方法によれば、ステップＳ７２以降、一つのバケットだけを対象とし、その中の実数から最近傍値を見つけることができる。このため、計算時間を大幅に短縮することができる。

上記実施形態では、前記最近傍値は、前記探索点ｑを超えずかつ該探索点に最も近い値のものが探索される。換言すれば、上記実施形態の方法では、前記最近傍値は、前記交点Ｃに最も近い物体モデル２内の格子点３ａが特定されることになる。

しかし、例えば、最近傍値は、前記探索点ｑを下回らない該探索点に最も近い値として定義されても良い。この場合、最近傍値は、前記交点Ｃに最も近い物体モデル２外の格子点３ａとして特定されることになる。この方法では、前記ラストフィルドバケットは、当該バケットよりも末尾バケット側で当該バケットに最も近いバケット番号を記憶させれば良い。また、ステップＳ７５については、「最大値」を「最小値」と読み替えれば良い。また、ステップＳ７３については、「最小値」を「最大値」と読み替える。さらに、ステップＳ７４については、x[ ipoint ] ≧ｑ＞x[ ipoint-1 ]となるx[ ipoint ] を最近傍値とすれば良い。

図１７、図１８には、本発明の他の実施形態の流体モデルが示される。
図１７の流体モデル３の格子点３ａは、ｒ−θ−ｚの円筒形空間のｒ方向（ｚ軸と直交する方向）をｓ分割、θ方向をｍ分割、ｚ方向をｎ分割したメッシュＭの交差点上に設けられている。また、前記直線Ｌは、前記メッシュに一致して定義される。即ち、
１）ｍ×ｎ個の直線がｒ方向と平行に配置、及び
２）ｓ×ｍ個の直線がｚ軸と平行に配置される。
これにより、流体モデル３の全ての格子点３ａが直線Ｌ上に位置するように定められる。

さらに、図１８の流体モデル３の格子点３ａは、ｒ−θ−φの球形空間のｒ方向をｓ分割、θ方向をｍ分割、φ方向をｎ分割したメッシュの交差点上に設けられ、前記直線は、前記メッシュに一致して定義される。即ち、ｍ×ｎ個の直線がｒ方向に配置される。これにより、流体モデル３の全ての格子点３ａが直線Ｌ上に位置するように定められる。

本発明の効果を確認するため、ゴルフボールモデル及び流体モデルを次の仕様で設定し、流体モデルの境界を計算する時間が測定された。コンピュータ装置には、ＨＰ社のワークステーションxw9300（搭載メモリ３２GB）が使用された。

＜ゴルフボールモデル＞
半径２１．３５mmを有し、表面は、三角形の要素（ラグランジュ要素）で形成されている。表面の要素の節点総数は３７５５５９８個である。

＜流体モデル＞
半径２２０mm、軸方向長さ４４０mmの円筒形空間とした。このサイズにより、液体モデルの最外側面の壁の影響を無くすことができる。メッシュの分割数（ｒ、θ、ｚ）は、次の水準とした。
流体モデル１：２５２×２５０２×１００
流体モデル２：５０２×３００２×１００

境界点のサーチアルゴリズムとして、実施例では、請求項１で示したように、流体モデルを貫く直線を用いた。また、比較のために、下記のプログラムを下記のウエブサイトからダウンロードし、ＡＮＮによる最近傍探索法との性能差を比較した。
ＡＮＮプログラム:
A Library for Approximate Nearest Neighbor Searching
David M. Mount and Sunil Arya Version 1.1.1 Release Date: Aug 4, 2006
URL http://www.cs.umd.edu/~mount/ANN/
計算時の設定：
validate off
stats query_stats
dim 3
data_size 150700400
query_size 3755598
read_data_pts test1-data.pts
read_query_pts test1-query.pts
bucket_size 1
near_neigh 1
split_rule suggest
shrink_rule none
build_ann
epsilon 0.0
run_queries priority
テストの結果は、表１に示される。なお、表１において、ＡＮＮ１は、バケットサイズをデフォルト（bucket_size 1）のまま使用したもの、ＡＮＮ３は、計算速度の向上のため、バケットサイズが３に設定されたものを示している。また、境界点の探索時間に加え、ソリッドマーキングに要した時間についても記載されている。実施例のソリッドマーキングは、直線上で認識された境界となる格子点の間に存在する格子点を一括して非流体領域としてマークした。他方、比較例１（比較例２は未実施）のものでは、公知のFlood fillアルゴリズムに従ってソリッドマーキングを行った。

テストの結果、本発明の実施例では、計算速度が大幅に向上していることが確認できる。

１ コンピュータ装置
２ 物体モデル
２ａ ゴルフボールモデル
３ 流体モデル
３ａ 格子点

Claims (8)

1. 三次元の固体を有限個の要素で表した物体モデルと、流体が流れる空間を規定する流体モデルとを用いて少なくとも前記流体の流れを解析する流体シミュレーションを行うに際し、
   ａ）コンピュータに、表面が複数の平面からなる前記物体モデルを入力するステップ、
   ｂ）前記コンピュータに、前記物体モデルを覆うように三次元空間に多数の流体解析用の格子点が定義された前記流体モデルを入力するステップ、
   ｃ）前記コンピュータに、前記流体モデルの前記格子点を通りかつ前記流体モデルを貫く複数の直線を定義するステップ、
   ｄ）前記コンピュータが、前記物体モデルの表面と前記直線との交点の座標値を計算しかつこの交点の座標値を記憶するステップ、及び
   ｅ）前記コンピュータが、上記ステップｄで得られた交点の座標値に基づいて、前記流体モデルの全ての格子点を固体領域又は流体領域のいずれかに決定するステップ
   を含んで前記流体モデルと前記物体モデルとの境界位置を前記流体モデルに定義し、
   前記ステップｅは、前記コンピュータの記憶部に、探索用データベースを作る前処理ステップと、
   前記コンピュータが前記探索用データベースを使用して前記交点の座標値に最も近い前記格子点の座標値の探索を行う探索ステップとを含み、
   前記探索用データベースは、前記直線上の格子点の座標値がなす実数群の区間（x[1]〜x[n]）を小領域であるバケットで隙間無くかつ等間隔でＱ個に区分されたバケット列を有し、
   前記各バケットには、先頭バケットから末尾バケットまで順番に割り当てられた固有のバケット番号と、
   各バケットの中に含まれる前記実数の個数であるバケットサイズとが記憶されるとともに、
   少なくとも前記先頭バケットを除いた各バケットには、前記バケットサイズがゼロでなくかつ当該バケットに最も近いバケットのバケット番号を示すラストフィルドバケットが記憶され、
   前記探索ステップは、前記交点の座標値が属するバケットを計算する第１のステップと、
   前記計算されたバケットのバケットサイズを参照する第２のステップと、
   前記参照されたバケットサイズが０でない場合、当該参照されたバケットから前記交点の座標値の最近傍値を計算する一方、前記参照されたバケットサイズが０の場合、そのバケットのラストフィルドバケットが示すバケットから前記交点の座標値の最近傍値を計算する第３のステップとを含むことを特徴とする流体シミュレーションにおける境界位置決定方法。
2. 前記流体モデルの前記格子点は、ｘ−ｙ−ｚの三次元空間のｘ軸方向をｓ分割、ｙ軸方向をｍ分割、ｚ軸方向をｎ分割したメッシュの交差点上に設けられ、
   前記直線は、前記メッシュに一致して定義される請求項１に記載の流体シミュレーションにおける境界位置決定方法。
3. 前記流体モデルの格子点は、ｒ−θ−ｚの円筒形空間のｒ方向をｓ分割、θ方向をｍ分割、ｚ方向をｎ分割したメッシュの交差点上に設けられ、
   前記直線は、前記メッシュに一致して定義される請求項１記載の流体シミュレーションにおける境界位置決定方法。
4. 前記流体モデルの格子点は、ｒ−θ−φの球形空間のｒ方向をｓ分割、θ方向をｍ分割、φ方向をｎ分割したメッシュの交差点上に設けられ、
   前記直線は、前記メッシュに一致して定義される請求項１記載の流体シミュレーションにおける境界位置決定方法。
5. 前記最近傍値は、前記交点の座標値に最も近い物体モデル内の格子点の座標値とし、
   前記探索用データベースには、
   前記ラストフィルドバケットとして、当該バケットよりも先頭バケット側で当該バケット番号に最も近いバケット番号が記憶されている請求項１記載の流体シミュレーションにおける境界位置決定方法。
6. 前記探索用データベースの少なくともバケットサイズが０でないバケットには、そのバケット内の実数のうちの最小値及び最大値がさらに記憶されている請求項１乃至５のいずれかに記載の流体シミュレーションにおける境界位置決定方法。
7. 前記第３のステップにおいて、参照されたバケットサイズが０でない場合、前記交点の座標値と、そのバケットの最小値とを比較し、前記交点の座標値が最小値よりも小さい場合には、そのバケットのラストフィルドバケットの最大値を最近傍値とする請求項６記載の流体シミュレーションにおける境界位置決定方法。
8. 前記第３のステップにおいて、参照されたバケットサイズが０でない場合、前記交点の座標値と、そのバケットの最小値とを比較し、前記交点の座標値が最小値以上の場合には、そのバケットに含まれる実数に対してx[ ipoint ] ≦ｑ＜x[ ipoint + 1 ]となるx[ ipoint ] を最近傍値とする請求項６又は７記載の流体シミュレーションにおける境界位置決定方法。