JP5544511B2 - 渦抽出装置、渦抽出方法、渦抽出プログラム、および渦抽出表示システム - Google Patents

Description

本件は、流体空間内の渦を抽出する渦抽出装置および渦抽出方法に関する。また本件は、プログラムを実行する演算装置により実行されてその演算装置を渦抽出装置として動作させる渦抽出プログラムに関する。さらに本件は、渦抽出装置とその渦抽出装置で抽出された渦を表示する渦表示装置とを有する渦抽出表示システムに関する。

流体空間内の渦を抽出して可視化する技術は、例えば、航空機の翼の周りの風の挙動、船のスクリューの周りの水流の挙動、人体の気管内を流れる空気の挙動等、様々な分野での流体の挙動を知る上で極めて重要である。

渦の可視化を実現するための渦の抽出の技術が従来より種々提案されている。この可視化技術における大きな課題は、渦の抽出を高精度に行なおうとすると厖大な量のデータを取り扱う必要があることである。このため、大量のデータを高速に取り扱うことのできる高速の大型演算機を必要とする。この渦抽出のための計算量を削減することも必要ではあるが、さらに大きな課題は、可視化の段階である。渦抽出の計算に用いたデータ量が膨大であることに対応して、抽出した渦を表わすデータのデータ量も厖大となり、渦抽出に用いた計算機と同レベルの高いスペックを持つ計算機を用いないと可視化処理が難しい。すなわち、スペック不足の計算機で可視化すると、表示された渦の回転、移動、拡大、縮小等の操作に対するレスポンスが遅く必要な表示がなかなか得られない。また、渦は、時々刻々と変化する。このため、渦の挙動を知るためには渦の時間的な変化を知る必要がある。しかしながらデータ量の厖大さに起因して、渦の時間的な変化の可視化が難しい。

高速な可視化処理を実現しようとしてデータを間引いて可視化すると、計算コストを投じて渦を正確に抽出したにもかかわらず、間引きにより重要な現象が失われて表示されずに見逃されてしまうおそれがある。

ここでは、渦抽出に関する文献のほか、厖大な量のデータを取り扱うための領域分割法の一例が示されている文献、および後述する説明で参考となる文献も挙げておく。

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

Ｄ．Ｓｕｊｕｄｉ， Ｒ．Ｈａｉｍｅｓ，ＩＤＥＮＴＩＦＩＣＡＴＩＯＮ ＯＦ ＳＷＩＲＬＩＮＧ ＦＬＯＷ ＩＮ ３−Ｄ ＶＥＣＴＯＲ ＦＩＥＬＤＳ，Ｐａｐｅｒ ｏｆ Ａｍｅｒｉｃａｎ Ｉｎｓｔｉｔｕｔｅ ｏｆ Ａｅｒｏｎａｕｔｉｃｓ ａｎｄ Ａｓｔｒｏｎａｕｔｉｃｓ，１９９５ Ｔ．Ｗｅｉｎｋａｕｆ，Ｃｏｒｅｓ ｏｆ Ｓｗｉｒｌｉｎｇ Ｐａｒｔｉｃｌｅ Ｍｏｔｉｏｎ ｉｎ Ｕｎｓｔｅａｄｙ Ｆｌｏｗｓ，Ｐａｐｅｒ ｏｆ ＩＥＥＥ，２００７ Ｏｎ ｔｈｅ ｉｄｅｎｔｉｆｉｃａｔｉｏｎ ｏｆ ａｖｏｒｔｅｘ， Ｊｉｎｈｅｅ Ｊｅｏｎｇ ａｎｄ Ｆａｚｌｅ Ｈｕｓｓａｉｎ，Ｊａｕｒｎａｌ ｏｆ ｆｌｕｉｄ ｍｅｃｈａｎｉｃｓ，ｖｏｌ．２８５，ｐｐ６４−９４，１９９５ Ｍａｒｃｈｉｎｇ ｃｕｂｅｓ：Ａ ｈｉｇｈ ｒｅｓｏｌｕｔｉｏｎ ３ Ｄ ｓｕｒｆａｃｅ ｃｏｎｓｔｒｕｃｔｉｏｎ ａｌｇｏｒｉｔｈｍ， Ｗｉｌｌｉａｍ Ｅ．Ｌｏｒｅｎｓｅｎ，Ｈａｒｖｅｙ Ｅ．Ｃｌｉｎｅ，Ｐｒｏｃｅｅｄｉｎｇｓ ｏｆ ｔｈｅ １４ｔｈ ａｎｎｕａｌ ｃｏｎｆｅｒｅｎｃｅ ｏｎ Ｃｏｍｐｕｔｅｒ ｇｒａｐｈｉｃｓ ａｎｄ ｉｎｔｅｒａｃｔｉｖｅ ｔｅｃｈｎｉｑｕｅｓ，Ｐａｇｅｓ：１６３−１６９，１９８７

本件開示の渦抽出装置等の課題は、重要な情報を保持し、かつ、通常のパーソナルコンピュータ等一般的なスペックの計算機等で高速な表示処理が可能な渦データを抽出することにある。

本件開示の渦抽出装置は、データ取得部と、第１の数値データ算出部と、テキストデータ作成部と、データ出力部とを有する。

データ取得部は、流体空間内に分散した複数の各点の流速データを取得する。

第１の数値データ算出部は、データ取得部で取得した流速データに基づいて、流体空間内を複数の第１のボクセルに分割したときの各第１のボクセルの各頂点ごとの、各頂点が流体の渦領域の内側にあるか外側にあるかを指標する第１の数値データを算出する。

テキストデータ作成部は、データ取得部で取得した流速データに基づいて流体の渦の中心線を求める。また、このテキストデータ作成部は、第１の数値データ算出部で算出された各頂点ごとの第１の数値データに基づいてその中心線の周りの渦領域の寸法を求める。このテキストデータ作成部は、このようにして求めた中心線と寸法とを表わすテキストデータを作成する。

データ出力部は、テキストデータ作成部で作成されたテキストデータを出力する。

本件開示の渦抽出プログラムおよび渦抽出方法は、本件の渦抽出装置の思想をプログラムおよび方法として把握したものである。

また本件開示の渦抽出表示システムは、本件の渦抽出装置と、その渦抽出装置から出力されたデータを受け取って渦を表示する渦表示装置を有する。

本件によれば、渦領域の中心線と寸法が相対的に高い分解能の流速データに基づいてテキストデータとして作成される。これらの重要な情報が先に作成されているために、その後、本件の渦抽出装置よりもスペックの低い装置で表示を行なうために、ボクセルデータやポリゴンデータ等の削減を仮に行なったとしても、重要な情報が失われることなく高速な渦表示処理が可能となる。

本件の渦抽出装置の一実施形態である大規模計算機を含む、システムの一例を示す模式図である。 本件の渦抽出プログラムの一実施形態を概括的に示した図である。 本件の渦抽出装置の一実施形態を概括的に示した図である。 本件の渦抽出方法の一実施形態の概要を示すフローチャートである。 大規模計算機の内部構成の概要を示す図である。 大規模計算機で実行される渦抽出プログラムの各ステップを示すフローチャートである。 内挿処理の説明図である。 ピクセル（ボクセル）の寸法の決め方のアルゴリズムの一例を示す図である。 領域分割処理の説明図である。 渦度の演算方法の説明図である 渦領域抽出処理の説明図である。 図６の渦の中心線分作成処理（Ｓ２８）の詳細フローを示すフローチャートである。 ボクセルの模式図である。 回転判定アルゴリズムの説明図である。 付加的に追加可能な回転判定アルゴリズムの説明図である。 ポリゴン数削減処理の説明図である。 渦表示の一例を示す図である。 渦を、中心線が左右に延びる方向から見て半透明で表示した図である。 中心線が紙面に垂直な方向に延びる方向から見て、渦を半透明で示した図である。 図１８の表示画面上にさらに中心線と渦の領域を表わす線を重畳して示した図である。 図１９の表示画面上にさらに中心線と渦の領域を表わす線を重畳して示した図である。 本実施形態の変形例を示した図である。

以下、本件の実施形態について説明する。

図１は、本件の渦抽出装置の一実施形態である大規模計算機を含む、システムの一例を示す模式図である。

ここでは、大規模計算機１０と、ファイルサーバ２０と、可視化用計算機３０が、インタネット等の通信回線網４０を介して接続されている。

大規模計算機１０は、本件の渦抽出装置の一実施形態に相当する計算機であり、大量のデータの高速演算処理が可能である。

ファイルサーバ２０は、大規模計算機１０での渦抽出の計算の基になる流速データや、大規模計算機１０での渦抽出の計算結果を表わすデータがファイリングされる装置である。

可視化用計算機３０は、渦表示装置の一例であって、大規模計算機１０での渦抽出の計算結果を表わすデータを受け取って表示画面３１上に表示する計算機である。この可視化用計算機３０は、一般に市販されているレベルのパーソナルコンピュータ（以下、「ＰＣ」と略記する）等の計算機であってもよい。

尚、この図１には、可視化用計算機３０を１台のみ示しているが、これは代表的に１台示したのであって、複数台あるいは多数台存在していてもかまわない。

ここでは、ファイルサーバ２０には、流体空間内の流速分布を表わす、その流体空間内に分散した複数の各点の流速データが格納されているものとする。

この流速データは、例えば風洞実験等を行なって実際に計測したデータであってもよく、あるいは、ここには図示しない、シミュレート用の計算機でシミュレートした結果得られたデータであってもよい。あるいは、図１に示す大規模計算機１０がシミュレート用の計算機を兼ねていてもよい。

ファイルサーバ２０内に格納されている流体空間内の流速分布を表わす流速データは、大規模計算機１０に入力され（矢印Ａ）、大規模計算機１０で渦抽出計算が行なわれる。その大規模計算機１０での渦抽出計算結果を表わすデータは、再度ファイルサーバ２０に送られてファイルサーバ２０に一旦格納される（矢印Ｂ）。その後、ファイルサーバ２０内の渦抽出計算結果を表わすデータは、可視化用計算機３０からの要求に応じて可視化用計算機３０に送られる（矢印Ｃ）。可視化用計算機３０ではその渦抽出計算結果を表わすデータを受け取ってそのデータに基づいて渦を表示する。

ファイルサーバ２０には、時々刻々と変化する流速分布のデータが時系列に蓄積されており、大規模計算機１０では各時刻ごとの渦の抽出が行なわれる。したがって可視化用計算機３０では、時々刻々と変化する渦の移り変わり（発生、消滅、合流、分岐等）の表示が可能である。また、ある一時点の渦を方向を変えて、あるいは拡大、縮小等を行なって表示することも可能である。

図２は、本件の渦抽出プログラムの一実施形態を概括的に示した図である。

この図２に示す渦抽出プログラム１００は、大規模計算機１０にインストールされてその大規模計算機１０で実行され、その大規模計算機１０を、本件の渦抽出装置の一例として動作させるプログラムである。この渦抽出プログラム１００は、データ取得部１１０、補間演算部１２０、第１の数値データ算出部１３０、第２の数値データ算出部１４０、テキストデータ作成部１５０、ポリゴン数削減処理部１６０、およびデータ出力部１７０で構成されている。これらの各部１１０〜１７０は渦抽出プログラム１００を構成する各プログラム部品である。各部１１０〜１７０の作用については後述する。

図３は、本件の渦抽出装置の一実施形態を概括的に示した図である。

この図３に示す渦抽出装置は、データ取得部２１０、補間演算部２２０、第１の数値データ算出部２３０、第２の数値データ算出部２４０、テキストデータ作成部２５０、ポリゴン数削減処理部２６０、およびデータ出力部２７０を有する。この図３に示す各部２１０〜２６０は、図２に示す渦抽出プログラム１００を構成する、図３と同じ名称を付した各プログラム部品１１０〜１７０が大規模計算機１０内で実行されることにより大規模計算機１０内に実現する機能をあらわしている。図２と図３とでは同じ名称を用いているが、図２の各部１１０〜１７０は、ソフトウェアのみを指し、図３の各部２１０〜２７０は、図２に示す同じ名称のプログラム部品と大規模計算機１０のハードウェアやＯＳ（Ｏｐｅｒａｔｉｎｇ Ｓｙｓｔｅｍ）等との複合により実現する機能を指している。したがって、ここでは、図３の各部２１０〜２７０の作用を説明することで、図２の各部１１０〜１７０の作用の説明を兼ねるものとする。

図３に示す渦抽出装置２００を構成するデータ取得部２１０は、図１に示すファイルサーバ２０から通信回線網４０を介して、渦抽出対象の流体空間内に分散した複数の各点の流速データを取得する機能を有する。

また、補間演算部２２０は、データ取得部２１０で取得した流速データに基づく補間演算により、流体空間を構造格子に分割したときの複数の第１のボクセルの各頂点ごとの流速データを算出する。

また、第１の数値データ算出部２３０は、流体空間内を複数の第１のボクセルに分割したときの各第１のボクセルの各頂点ごとの、各頂点が流体の渦領域の内側にあるか外側にあるかを指標する第１の数値データを算出する。この第１の数値データ算出部２３０は、この第１の数値データを、データ取得部２１０で取得した流速データに基づいて算出する。

さらに、第２の数値データ算出部２４０は、少なくとも、第１の数値データ算出部２３０で算出された第１の数値データに基づいて判定される渦領域を構成する各第１のボクセルそれぞれについて、各第１のボクセルごとの渦の強さを指標する第２の数値データを算出する。この第２の数値データ算出部２４０は、この第２の数値データを、データ取得部２１０で取得した流速データに基づいて算出する。この第２の数値データ算出部２４０では、各ボクセルごとの第２の数値データを、渦領域の内部に限って算出してもよく、渦領域の内部であるか外部であるかの判定なしに、流体空間全域に亘って算出してもよい。

また、テキストデータ作成部２５０は、データ取得部２１０で取得した流速データに基づいて流体の渦の中心線を求める。また、このテキストデータ作成部２５０は、第１の数値データ算出部２３０で算出された各頂点ごとの第１の数値データに基づいてその中心線の周りの渦領域の寸法を求める。そしてこのテキストデータ算出部２５０は、求めた中心線と寸法とを表わすテキストデータを作成する。

さらに、ポリゴン数削減処理部２６０は、第１のデータ算出部２３０で算出された第１の数値データの等値面で構成される、渦領域の表面を画定するポリゴンの数を計数する。そして、複数の第１のボクセルを、ポリゴンの数が設定値（図３に示す「設定値」２６１）以下となるまで、流体空間が粗く分割された第２のボクセルに統合する。すなわち、第１の数値データを、第２のボクセルの各頂点に対応する第１の数値データにまで間引く。このポリゴン数削減処理部２６０はさらに、第２の数値データ算出部２４０で算出された第２の数値データに基づいて、統合された第２のボクセルそれぞれに対応する第２の数値データを作成する。ここで、「設定値」２６１は、可視化用計算機３０（図１参照）からあらかじめ、その可視化用計算機３０ごとに個別に指定された値であってもよい。あるいはこの「設定値」２６１は、平均的なＰＣで高速な可視化処理が可能なように大規模計算機１０側で、可視化用計算機３０の種類等に寄らずに一律に設定した値であってもよい。あるいはこの「設定値」２６１は、図２に示す渦抽出プログラム１００が定数としてもっている値であってもよい。

さらに、データ出力部２７０は、テキストデータ作成部２５０で作成されたテキストデータを出力する。またこのデータ出力部２７０は、ポリゴン数削減処理部２６０で得られた第２のボクセルの各頂点に対応する第１の数値データも出力する。このデータ出力部２７０は、テキストデータおよび第１の数値データに加え、さらに、第２のボクセルに対応する第２の数値データも出力する。

前述した通り、この渦抽出装置２００から出力されたデータはファイルサーバ２０に一旦格納され、その後可視化用計算機３０に送られる。可視化用計算機３０では、それらのデータに基づいて表示画面３１上に渦が表示出力される。

図４は、本件の渦抽出方法の一実施形態の概要を示すフローチャートである。

この図４に示す渦抽出方法は、図１に示す大規模計算機１０内での図２に示す渦抽出プログラム１００の実行により実施される渦抽出方法である。この渦抽出方法は、データ取得過程（Ｓ１１）、補間演算過程（Ｓ１２）、第１の数値データ算出過程（Ｓ１３）、第２の数値データ算出過程（Ｓ１４）を有する。この渦抽出方法は、さらに、テキストデータ作成過程（Ｓ１５）、ポリゴン数削減処理過程（Ｓ１６）、データ出力過程（Ｓ１７）を有する。

これらの各過程Ｓ１１〜Ｓ１７における実行内容は、図３に示す渦抽出装置２００の各部２１０〜２７０の作用とそれぞれ同一であり、ここでは重複説明は省略する。

次に更に詳細な実施形態について説明する。以下に説明する実施形態においてもシステムの全体構成は、図１をそのまま流用する。

図５は、図１に示す大規模計算機の内部構成の概要を示す図である。

この図５に示す大規模計算機１０は、プログラムを実行する複数（ここでは４台）のＣＰＵ（Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）１１Ａ〜１１Ｄと、それらのＣＰＵ１１Ａ〜１１Ｄで共用されるメモリ１２を有する。これら複数のＣＰＵ１１Ａ〜１１Ｄは、流体空間を分割した各領域における渦の抽出を担当している。ここでは、このように複数のＣＰＵで流体空間内の領域を分けて渦抽出処理を行なっており、高速の演算処理を可能としている。

図６は、図１，図５に示す大規模計算機１０で実行される渦抽出プログラムの各ステップを示すフローチャートである。

ここでは、図６に示すデータ入力処理（Ｓ２１）からデータ出力処理（Ｓ２９）までの各処理が実行される。

ここで、一例としてデータ入力処理（Ｓ２１）、内挿処理（Ｓ２２）、領域分割処理（Ｓ２３）、およびデータ出力処理（Ｓ３０）は、図５に示す複数のＣＰＵ１１Ａ〜１１Ｄのうちの１つのＣＰＵ１１Ａが担当している。その他の処理（Ｓ２４〜Ｓ２９）は、他の複数のＣＰＵ１１Ｂ〜１１Ｄが、流体空間を複数に分けたときの各領域についての計算を分担している。図１，図５に示す大規模計算機１０には、時々刻々変化する時系列データとしての流速データが入力される。ある時刻についての流速データに関するＣＰＵ１１Ａでの領域分割処理（Ｓ２３）が終わると、ＣＰＵ１１Ｂ〜１１Ｄではそれぞれ分担した領域についてその時刻についての渦抽出処理が実行される。その間、ＣＰＵ１１Ａでは、ＣＰＵ１１Ｂ〜１１Ｄが分担した領域について渦抽出処理を実行している間に次の時刻の流速データのデータ入力処理（Ｓ２１）等が行なわれる。ＣＰＵ１１Ｂ〜１１Ｄである時刻についての渦抽出処理（Ｓ２３〜Ｓ２９）が終了すると、ＣＰＵ１１Ａがその抽出処理結果を出力する（Ｓ３０）。このように、この大規模計算機１０では、複数のＣＰＵ１１Ａ〜１１Ｄで時々刻々変化する流速データについて役割りを分担して並行処理が行なわれる。

図６に示すフローチャート中のデータ入力処理（Ｓ２１）は、前述の実施形態（図３に示す渦抽出装置２００で代表させる）のデータ取得部２１０に対応する処理を実行する。すなわち、このデータ入力処理（Ｓ２１）は、大規模計算機１０にファイルサーバ２０に格納されている時系列データとしての流速データを順次受信させ、その受信した流速データをこの図６に示す渦抽出プログラムに取り込む。

また、内挿処理（Ｓ２２）は、大規模計算機１０に、図３に示す渦抽出装置２００の補間演算部２２０に対応する処理を実行させる。

図７は、内挿処理の説明図である。

図７（Ａ）は、内挿処理前の非構造格子の一例を示した図、図７（Ｂ）は、図７（Ａ）の非構造格子に構造格子を重ねて示した図、図７（Ｃ）は内挿処理後の構造格子を示した図である。これら図７（Ａ）〜（Ｃ）それぞれの中央には、飛行機の翼を模擬したときのその翼の断面が示されている。尚、流体空間（ここでは翼の周囲の空間）は三次元であるが、ここでは簡単のため、二次元的に示してある。

図７（Ａ）には、不定形の多角形からなるモザイク模様が示されている。このモザイク模様を構成する多角形は、周辺ほど平均的な面積が広く（平均的な辺の長さが長く）、中央の翼に近づくにつれて平均的な面積が狭く（平均的な辺の長さが短かく）なっている。このモザイク模様を構成する各多角形の各頂点には、この流体空間内を流れる流体（例えば空気）の、各頂点の位置における流速を表わす流速データが対応づけられている。この流速データは、その対応づけられている位置の、３次元的な流速方向と大きさからなる流速ベクトルを表わしている。

この流速データは、この流体空間内の分散した複数の各点の流速を表わしてはいるが、このままでは取扱いが不便である。そこで、図７（Ｂ）に重ねて示すように、ｘｙの２方向（実際はｘｙｚの３方向）に整然と区画された４角形（６面体）のピクセル（ボクセル）を考える。そして、図７（Ａ）に示す不定形の多角形の各頂点に対応づけられた流速データに基づく内挿処理（補間演算処理）により、各ピクセル（各ボクセル）の各頂点に対応づけられた流速データを求める。こうすることにより、図７（Ｃ）に示すように、流体空間内に等間隔に配列された各点の流速データが求められて、その後の演算処理が容易となる。

各ピクセル（各ボクセル）の１辺の長さ（流体空間内での流速データの間隔）は、特に限定されるものではないが、データ間隔が粗過ぎると重要な現象が抽出されないおそれがある。一方、データ間隔が密過ぎると、大規模計算機１０の処理能力が高くても演算に多大な時間がかかるおそれがあり、またいくら密であってももともとの非構造格子（図７（Ａ））の流速の分解能以上の分解能は得られない。したがって非構造格子（図７（Ａ））を構成する多数の多角形の平均的な面積（３次元では平均的な体積）と同程度の面積（体積）のピクセル（ボクセル）となるように、ピクセル（ボクセル）の寸法を決めることが好ましい。あるいは、多数の多角形の平均的な辺の長さと同程度の長さの辺を持つピクセル（ボクセル）となるようにピクセル（ボクセル）の寸法を決めてもよい。

図８は、ピクセル（ボクセル）の寸法の決め方のアルゴリズムの一例を示す図である。

図８の横軸は、図７（Ａ）の非構造格子に出現する多角形の辺の長さの最小値ｍｉｎと最大値ｍａｘとの間を、対数軸上で１０等分したものである。また図８の縦軸は、横軸（辺の長さ）を対数軸上で１０等分したときの各区画ごとの辺の出現数である。

図７（Ａ）に示す非構造格子を構成する多角形の辺の長さを調べて図８の横軸の各区画ごとの辺の長さの出現数のうち、最小の区画Ａと最大の区画Ｂの各出現数を計数する。最小の区画Ａの出現数をａ、最大の区画Ｂの出現数をｂとし、最小の区画Ａの辺の長さの代表値をＬｍｉｎ、最大の区画Ｂの辺の長さの代表値をＬｍａｘ、図７（Ｃ）に示す構造格子（ボクセル）の一辺の長さをＬとしたとき、（１）式

に従って構造格子（ボクセル）の１辺の長さＬを算出してもよい。例えばこのようにして算出した１辺の長さＬのボクセルの各頂点の流速データが内挿処理（補間演算処理）により求められる。ただし、ここに示した長さＬの求め方は一例であり、非構造格子の辺の長さの単純平均、流体空間の一部のサンプル領域の辺の平均的な長さ等、様々な算出方法のいずれを採用してもよい。例えば、図７（Ａ）の場合、中央の翼に近づくほど平均的な面積が小さい（辺の平均的な長さが短い）非構造格子となっている。そこで、翼の近傍に関心があるときは、翼の近傍領域の非構造格子の平均的な辺の長さを求めて、その平均的な長さの辺を持つ構造格子となるように内挿処理（補間演算処理）を行なってもよい。

尚、内挿処理（補間演算処理）自体は良く知られた手法であり、ここでは、内挿処理（補間演算処理）自体についての説明は省略する。また、ここでは非構造格子としての流速データが入力されるものとして説明したが、もともとの流速データが図７（Ｃ）に示すような構造格子の流速データであったときは、この内挿処理（Ｓ２２）は不必要であるためスキップされる。

上記の内挿処理（Ｓ２２）が終了すると、領域分割処理（Ｓ２３）を経て次に流速勾配算出処理（Ｓ２３）が実行される。前述した通り、図６のデータ入力処理（Ｓ２１）および内挿処理（Ｓ２２）は、図５に示す複数のＣＰＵ１１Ａ〜１１Ｄのうちの１つのＣＰＵ１１Ａで実行される。これに対し、流速勾配算出処理（Ｓ２４）は、他の複数のＣＰＵ１１Ｂ〜１１Ｄで、流体空間が複数に分割された各領域について分担して実行される。そのため、ここでは領域分割処理（Ｓ２３）が行なわれる。

図９は、領域分割の説明図である。ここでも、簡単のため、２次元的な構造格子で示してある。

ここでは、流体空間Ｄが、一部が重なるようにして３つの領域Ｄ１〜Ｄ３に分割されている。斜線を付した２つの領域Ｄ１２，Ｄ２３のうちの領域Ｄ１２は、領域Ｄ１と領域Ｄ２の双方に含まれており、もう１つの領域Ｄ２３は、領域Ｄ２と領域Ｄ３の双方に含まれている。このように複数の領域に含まれる領域を袖領域と称する。

図５に示す複数（３つ）のＣＰＵ１１Ｂ〜１１Ｄは、各領域Ｄ１〜Ｄ３について、図６のフローチャートにおける局所固有値算出処理（Ｓ２４）以降の渦抽出処理を実行する。ここで、図９に示すように袖領域Ｄ１２，Ｄ２３が設定されているため、袖領域Ｄ１２，Ｄ２３については、その袖領域Ｄ１２，Ｄ２３を含む領域Ｄ１〜Ｄ３の処理を分担している複数のＣＰＵでほぼ同じ処理結果を得ることができる。したがって、処理終了後に流体空間Ｄの全域の処理結果として円滑に統合することができる。

以下では、分割された領域であるか、流体空間全域であるかを、特に必要なとき以外は区別せずに説明を続行する。

図６のフローチャートでは、流速勾配算出処理（Ｓ２４）、局所固有値算出処理（Ｓ２５）、および渦領域抽出処理（Ｓ２７）の組合せが、図３の第１の数値データ算出部２３０での処理に対応する。ただし、渦領域抽出処理（Ｓ２７）は、渦領域の抽出制度を向上させるための処理であり、この渦領域抽出処理（Ｓ２７）を必ずしも行なわなくても第１の数値データ算出処理としては十分である。

流速勾配算出処理（Ｓ２４）では、以下の（２）式に基づいて、各ボクセルの各頂点の流速勾配が算出される。

ここでは、流体空間内のｘ，ｙ，ｚの各方向の速度をそれぞれｕ，ｖ，ｗで表わし，ｕ，ｖ，ｗを代表させてｑで表わし、ｘ，ｙ，ｚを代表させてｒで表わす。また、各頂点を代表させて、その頂点の座標を（ｉ，ｊ，ｋ）で表わす。さらに、ｉ，ｊ，ｋを代表させてｎで表わし、そのｉ，ｊ，ｋやｎを省略することがある。

このとき、頂点（ｉ，ｊ，ｋ）の流速勾配は、

であらわされる。ただし、Δｒはｒ方向のボクセルの１辺の長さを表わしている。

ここで、図７の中央に示す翼との境界や流体空間を画定する壁面等、停止した物体との境界上ではｑ_ｒ＝０となるように処理を行なうことが好ましい。

ここでは、（２）式１つで示しているが、この（２）式は、下記の（３）式の流速行列Ｊの９つの各要素を表わしている。

各ボクセルの各頂点（ｉ，ｊ，ｋ）について、上記の（２）式による流速勾配の算出が行なわれると、次に局所固有値算出処理（Ｓ２５）が行なわれる。

この局所固有値算出処理（Ｓ２５）では、上記の（２）式に従って求めた流速勾配を要素とする流速行列Ｊを、

とし、

ただし、Ｊ^ＴはＪの転置行列である。
としたときの行列Ａ

の３つの固有値λ_ａ，λ_ｂ，λ_ｃを算出する。

これらの３つの固有値λ_ａ，λ_ｂ，λ_ｃが算出されると、これら３つの固有値λ_ａ，λ_ｂ，λ_ｃをその値順に並べ、λ_１＜λ_２＜λ_３とする。ただし、λ_１，λ_２，λ_３は、それぞれλ_ａ，λ_ｂ，λ_ｃのいずれかである。

ここでは、このようにして値順に並べた３つの固有値λ_１，λ_２，λ_３のうちの中央の値を持つ固有値λ_２を抽出する。このλ_２は、図３の第１の数値データ算出部２３０で算出される第１の数値データの一例に相当する。

この固有値λ_２を各ボクセルの各頂点（ｉ，ｊ，ｋ）のそれぞれについて算出する。ある頂点（ｉ，ｊ，ｋ）について算出した固有値λ_２がλ_２＜０のときは、その頂点（ｉ，ｊ，ｋ）は渦領域の内部であると判定され、λ_２≧０のときはその頂点（ｉ，ｊ，ｋ）は渦領域の外部であると判定される。

この渦領域の判定アルゴリズムは、渦を大きく捉えたとき、渦は一次元的に延びる中心線の周りを周回する二次元的な動きであることに着目した判定アルゴリズムである。このアルゴリズム自体は、既知のアルゴリズムであり、Ｊｅｏｎｇらのλ_２法と称される（前掲の非特許文献３参照）。

ここでは、ある頂点（ｉ，ｊ，ｋ）についてλ_２＜０のが算出されたときは、そのλ_２をその値のままその頂点（ｉ，ｊ，ｋ）に対応づける。一方、ある頂点（ｉ，ｊ，ｋ）についてλ_２≧０のλ_２が算出されたときは、λ_２をλ_２＝０に置き換えた上で、そのλ_２＝０をその頂点に対応づける。

これによりλ_２＜０領域が渦領域、λ_２＝０の領域は渦の外の領域であるという、渦領域の大まかな抽出が行なわれる。

次に渦度絶対値算出処理（Ｓ２６）が実行される。この渦度絶対値算出処理（Ｓ２６）では渦度ωの絶対値｜ω｜が算出される。この渦度絶対値算出処理（Ｓ２６）は、図３の第２の数値データ算出部２４０での処理に対応する。

ここでは、この渦度ωは、ボクセルごとに算出される。

この渦度ωは、あるボクセルの流速ベクトルをＵとしたとき、

但し、∇は、演算子

であり、×は外積（ベクトル積）を表わしている。

で表記すると、

となる。

ここで、Ｕはボクセルを代表する流速ベクトルであるところ、流速データはボクセルの各頂点の流速データである。１つの手法として各頂点の流速データを基に内挿処理（補間演算処理）によりボクセルの中央の流速データを算出し、その算出した流速データをＵとして（９）式の演算を行なうことが考えられる。ただしここでは、内挿演算と、（９）式によるωの演算とを分けることなく、以下のようにしてωを算出する。ここでは、理解の容易さのため、２次元で説明する。

（９）式を２次元で書き直すと、

となる。

図１０は、渦度ωの演算方法の説明図である。ここでは、簡単のためｘ−ｙの２次元平面について示してある。

あるボクセル（ここでは２次元のため、あるピクセル）各頂点（ｉ，ｊ，ｋ），（ｉ＋１，ｊ，ｋ），（ｉ，ｊ＋１，ｋ），（ｉ＋１，ｊ＋１，ｋ）のｘ方向、ｙ方向の速度を、図１０に示すように、それぞれｕ，ｖにその頂点の座標（例えばｉ，ｊ，ｋ等）を付して表記する。

このとき、このピクセルの中央の点の、（１０）式に従う、渦度

は、

となる。実際には、（１０），（１１）式に従う２次元の演算ではなく、（９）式に従う３次元の演算が行なわれ、各ボクセルごとの渦度ωが算出され、その絶対値｜ω｜で表わされる渦の強さが求められる。この渦の強さ｜ω｜は、図３の第２の数値データ算出部２４０で算出される第２の数値データの一例に相当する。

次に、本実施形態では、図６のフローチャートの渦領域抽出処理（Ｓ２７）により、局所固有値算出処理（Ｓ２５）で一旦求めた、渦領域の内部（λ_２＜０）であるか外部（λ_２＝０）であるかを示す固有値λ_２が必要に応じて修正される。

図１１は、渦領域抽出処理の説明図である。

図１１（Ａ）は、λ_２＜０の領域を模式的に示した図、図１１（Ｂ）は、｜ω｜の値による渦領域分離の説明図である。ここでも簡単のため２次元平面で示してある。

図１１（Ａ）では、斜線を施した領域がλ_２＜０の領域（渦領域内部）であるとする。ここには、この段階では未だ不明ではあるが、図示の２つの渦が存在しているものとする。ここでは、局所固有値算出処理（Ｓ２５）で算出された固有値λ_２がλ_２＜０の領域内について、渦度絶対値算出処理（Ｓ２６）で求められた渦度絶対値、すなわち渦の強さ｜ω｜の値を調べる。そして｜ω｜がある閾値と比較され、閾値以下の｜ω｜を持つボクセルのλ_２をλ_２＝０に置き換える。例えばここでは図１１（Ｂ）に斜線を付して示した領域内のボクセルの｜ω｜が閾値以下であるとする。このように、λ_２を用いて一旦抽出した渦領域を｜ω｜を用いて必要に応じて分離することにより、渦領域がより正確に抽出される。

このようにして渦領域の抽出を行なった後、次に図６の渦の中心を示す線分作成処理（Ｓ２８）が実行される。この渦の中心を示す線分作成処理（Ｓ２８）は、図３のテキストデータ作成部２５０での処理に対応する。

図１２は、図６の渦の中心を示す線分作成処理（Ｓ２８）の詳細フローを示すフローチャートである。また、図１３は、ボクセルの模式図である。

図１３に示すように、ボクセルの各頂点には、図１３に矢印で示すようなベクトルで表現される流速データが対応づけられている。ここでは代表的に、１つの頂点の流速データを（ｕ，ｖ，ｗ）（ｕはｘ方向の流速、ｖはｙ方向の流速、ｗはｚ方向の流速）であらわす。

図１２に示す渦の中心を示す線分作成処理では、このボクセルの、ｚ方向を向いた２つの面（Ｉ面とＩＩ面）、ｙ方向を向いた２つの面、およびｘ方向を向いた２つの面の順に、以下の処理が実行される。ここでは、ｚ方向を向いた２つの面（Ｉ面とＩＩ面）について説明する。ｙ方向、ｘ方向を向いた各２つの面についても、方向が異なるだけであって、同じ処理が実行される。

ｚ方向を向いた２つの面（Ｉ面とＩＩ面）のうちの、先ずはＩ面について、そのＩ面の４つの隅に相当する４つの頂点の各流速ベクトルをＩ面に写像し、その写像したベクトルを取り扱う。換言すると、４つの各頂点の流速ベクトル（ｕ，ｖ，ｗ）のうち、ｚ方向分ｗを無視し、ｘ，ｙの成分（ｕ，ｖ）に基づいて流体が回転しているか否かが判定される。回転しているか否かの判定基準については後述する。Ｉ面について流体が回転していることが判定されると、次に同じｚ方向を向いたもう１つの面であるＩＩ面について同様の判定が行なわれる。すなわち、そのＩＩ面の４つの隅に相当する４つの頂点の各流速ベクトル（ｕ，ｖ，ｗ）のｚ成分２は無視され、ｘ，ｙ成分（ｕ，ｖ）に基づいて、流体が、Ｉ面での回転方向と同じ方向に回転しているか否かが判定される。さらに、後述する湧出し、吸込みについてもＩ面とＩＩ面とで同一か否かが判定される。Ｉ面とＩＩ面が共に同じ方向に回転しており、湧出し、吸込みについても同一であると判定されると、このボクセルに、正の定数Ｖｃを割り当てる。この定数Ｖｃは、このボクセルが渦の中心軸上のボクセルであることをあらわしている。これと同じ処理を、流体空間内の、少なくとも１つの頂点がλ_２＜０（渦領域内）であるボクセルについて実行する。さらにこれと同じ処理を、ｙ方向を向いた２つの面について実行し、さらに同じ処理をｘ方向を向いた２つの面について実行する。

図１２は、基本的には以上の処理をあらわしたフローチャートであり、このフローチャートに沿って処理内容を説明すると以下の通りとなる。

ステップＳ３１では，ｚ，ｙ，ｘのいずれの方向について処理を実行するか設定する。上述の説明はｚ方向についての処理である。

ステップＳ３２では、ｚ，ｙ，ｘの３方向とも処理が終了しているか否かが判定され、ｚ，ｙ，ｘの３方向とも処理が終了すると、ステップＳ４２に進む。ここでは、ステップＳ４２の処理の説明は後回しにする。ここでは、ステップＳ３１でｚ方向が設定されたものとして以下説明する。ステップＳ３３では、未処理のボクセルの中から１つのボクセルを次の処理実行対象のボクセルとして設定する。ステップＳ３４では、ステップＳ３１で設定した方向（ここではｚ方向）に関し、全てのボクセルの処理が終了したか否かが判定される。このステップＳ３４で現在設定されている方向（ここではｚ方向）に関し全てのボクセルの処理が終了したと判定されるとステップＳ３１に戻り、次の方向（例えばｚ方向の次はｙ方向）が設定される。ステップＳ３３で未処理のボクセルが設定されるとステップＳ３４を経由してステップＳ３５に進む。このステップＳ３５では、ステップＳ３３で設定されたボクセルの８つの頂点のうちのいずれか１つの頂点でもλ_２＜０を満足するか否かが判定される。８つの頂点の全てがλ_２＝０（前述の通り、λ_２≧０はλ_２＝０に置き換えてある）のときは、そのボクセルは渦領域の外にあるため、そのボクセルについては処理は不要であり、ステップＳ３３に戻って次のボクセルを設定する。

ステップＳ３５で８つの頂点のうちの１つでもλ_２＜０であると判定されると、ステップＳ３６に進み、現在設定されている方向（ここではｚ方向）を向いた２つの面（Ｉ面とＩＩ面）のうちの、Ｉ面について、流体が回転をしているか否かが判定される。この判定に当たっては、上述の通り、そのＩ面の４つの隅に相当する４つの頂点の流速ベクトルの、現在設定されている方向（ここではｚ方向）の流速成分は無視され、残りの２方向の流速成分が考慮される。回転の判定基準の説明はさらに後に回す。

Ｉ面について回転していないと判定された場合は、そのボクセルのその方向（ここではｚ方向）についての処理を中止してステップＳ３３に戻る（ステップＳ３７）。Ｉ面について回転していると判定されると、今度は同じ方向（ここではｚ方向）を向いたもう１つの面（ＩＩ面）について同様の判定が行なわれる（ステップＳ３８）。ＩＩ面について回転していないと判定された場合（ステップＳ３９）、およびＩＩ面について回転しているがＩ面についての回転方向と同一方向（後述する湧出し、吸込みについても同一）の回転ではないと判定された場合（ステップＳ４０）は、ステップＳ３３に戻る。ＩＩ面がＩ面と同じ方向（湧出し、吸込みについても同一）に回転していると判定されると、ステップＳ４１に進み、現在の処理対象のボクセルに、中心軸上のボクセルであることを表わす値Ｖｃを設定した上でステップＳ３３に戻る。ステップＳ３３では、現在処理を行なっている方向（ここではｚ方向）について未処理のボクセルの中から次の１つのボクセルを次に処理すべきボクセルとして設定する。

以上の処理がｚ，ｙ，ｘの３方向の全てについて終了すると、ステップＳ４２に進む。ステップＳ４２では、中心軸上のボクセルであることを示す値Ｖｃが設定されたボクセルのうち、ほぼ１次元的に並んだボクセル群にスプライン関数を当て嵌めて渦の中心線を算出する。ほぼ１次元的に並んだボクセル群という条件を入れているのは、１次元的に並んだボクセル以外に、それらのボクセルから離れた位置に別の副次的な中心線が存在することもあり得るからである。渦の中心線が求められると、次にステップＳ４３に進む。ここでは、中心線の長手方向中央の点で、その中心線と垂直に交わる平面を考えたとき、その平面上における渦領域（λ_２＜０の領域）の表面までの、中心線（と平面との交点）からの最短距離ｒ_ｍｉｎと最長距離ｒ_ｍａｘが求められる。これらの最短距離ｒ_ｍｉｎと最長距離ｒ_ｍａｘとにより渦領域を概括的に示した領域が定まり、その領域が中心線に沿って延びた渦の中心を示す線分が作成される。

図６の渦の中心を示す線分作成処理（Ｓ２８）では、以上のようにして求めた渦の中心線と渦の領域とによる渦の中心を示す線分のデータが作成される。この渦の中心を示す線分のデータは、中心線を表わす関数と、ｒ_ｍｉｎ，ｒ_ｍａｘからなる渦領域を表わすデータからなり、個々のポリゴンには対応づけられていないテキストデータとして作成される。

次に、図１２のステップＳ３６，Ｓ３８における、回転しているか否かの判定処理について説明する。

図１４は、回転判定アルゴリズムの説明図である。ここでも簡単のため、各ボクセルを２次元平面上に写像したときのピクセルとして示してある。

図１４（Ａ），（Ｂ）とも、斜線が付されたボクセルが流体が回転していると判定されるボクセルであり、その右隣りのボクセルは流体が回転していないと判定されるボクセルである。図１４（Ａ）では、斜線が施されたピクセルの４つの隅の流速ベクトルは、そのピクセルの各辺との角度９０°以内で各辺の外側を各辺に沿って同じ回転方向（ここでは時計回りの回転方向）を向いている。この場合、このピクセルの中に渦の湧出しの中心点が存在することになる。

また図１４（Ｂ）では、斜線が施されたピクセルの４つの隅の流速ベクトルは、各辺の内側を向き、各辺から９０°以内で各辺に沿って同じ回転方向を向いている。この場合このピクセルの中に渦の吸込みの中心点が存在することになる。

尚、図１４（Ａ），（Ｂ）の斜線が付されたピクセルの右隣りのピクセルは、上記の湧出し、吸込みのいずれの条件にも合致しない。

ここでは、同一方向（例えばｚ方向）を向いたＩ面とＩＩ面のうちの一方の面についてのみ説明したが、同一方向を向くＩ面とＩＩ面の双方の面が、同一方向に、湧出し吸込みも同一に回転していることをもって、そのピクセル（ボクセル）に、渦の中心線上のボクセルであることを表わす値Ｖｃが設定される。

但し、図１４を参照して説明した回転判定アルゴリズムでは、中心線上のボクセルであるにもかかわらず中心線上のボクセルであると判定されないことがまれに発生する。しかしながら、ここではスプライン関数をあてはめて中心線を定めているので、中心線上のボクセルにまれに値Ｖｃが設定されていなくてもスプライン関数が連続するため問題はない。

但し、図１４を参照して説明した判定アルゴリズムに、さらに以下の判定アルゴリズムを追加することで、中心線上のボクセルであることの判定精度を高めることも可能である。

図１５は、付加的に追加可能な回転判定アルゴリズムの説明図である。

図１５（Ａ），（Ｂ）とも、右下のピクセル（ボクセル）が流体が回転しているピクセル（ボクセル）であると判定される。ただしここではＩ面とＩＩ面のうちの一方の面のみ示してあり、Ｉ面とＩＩ面との双方で回転の条件を満たす必要があることは上述の通りである。

図１５（Ａ），（Ｂ）の右下の４角形のピクセルは、いずれも、図１４（Ａ），（Ｂ）のいずれの回転判定条件を満足していない。そこで、そのピクセルの対角を結んで三角形を作る。一対の対角を結ぶことで２つの三角形が形成され、別の一対の対角を結ぶことでもう２つの三角形が形成される。それら４つの三角形のそれぞれについて、図１４を参照して説明した判定基準と同じ判定基準を適用する。すなわち、図１５（Ａ），（Ｂ）のそれぞれに斜線を付した三角形に着目すると、その三角形の３つの頂点の流速ベクトルは各辺から９０°以内で各辺に沿って同一の回転方向を向いており、かつ湧出し、吸込み（図１５（Ａ），（Ｂ）はいずれも湧出し）の条件も満たしている。このような三角形が存在するボクセルにも、値Ｖｃを設定する。図１４の判定アルゴリズムに、この図１５の判定アルゴリズムを追加することで、渦の中心線上のボクセルであることの判定精度を高めてもよい。

次に、図６のフローチャート中の、ポリゴン数削減処理（Ｓ２９）について説明する。

図１６は、ポリゴン数削減処理の説明図である。ここでも、簡単のため、ポリゴンを２次元的なピクセルとして示してある。またここでは、これも簡単のため、縦４個×横４個の合計１６個のポリゴン（ピクセル）のみ示してある。

各ポリゴン（ピクセル）の各頂点には、λ_２が対応づけられている。ここでは、このλ_２の添字の２は省略し、座標を表わすｉ，ｊ等を添字として付してある。ただし、各座標のλ_２を代表させるときは添字２を付してλ_２と表記する。前述の通り、各λ_２は、λ_２＜０又はλ_２＝０である。ここでは、絶対数の小さいマイナスの値、例えば−０．１等を設定し、Ｍａｒｃｈｉｎｇ Ｃｕｂｅ法を適用してλ_２＝−０．１の等値面からなるポリゴンを生成する。このポリゴンは、図１６上では、線分としてあらわれている。Ｍａｒｃｈｉｎｇ Ｃｕｂｅ法自体は、例えば前掲の非特許文献４に記載されている公知の技術である。このλ_２＝−０．１の等値面からなるポリゴンは、渦領域の表面を画定するポリゴンである。図１６では、斜めの実線の線分が各ポリゴンに相当する。３次元ボクセルが配列されている流体空間では、１つずつのポリゴンは三角形となる。

このようにして生成したポリゴンの数を計数する。ここでは、図５を参照して説明したように複数のＣＰＵ１１Ｂ〜１１Ｄで領域を分担して（図９参照）処理を実行している。このためここでは、ＣＰＵ１１Ａが、それら複数のＣＰＵ１１Ｂ〜１１Ｄからそれぞれの分担領域のＤ１〜Ｄ３ポリゴン数の報告を受け、袖領域Ｄ１２，Ｄ２３におけるポリゴン数が二重に計数されないようにしながらポリゴン数を集計する。そしてこのポリゴン数が設定値を上回わった時は、各ＣＰＵ１１Ｂ〜１１Ｄに指示してボクセルを、ポリゴン数が設定値以下となるまで、流体空間が粗く分割されて新たなボクセルに統合させる。具体的には、例えば図１６に示す縦４個×横４個のボクセルを隣接した縦２個×横２個ボクセルを統合して、図１６に一点鎖線で示す新たなボクセルとする。この新たなボクセルに統合するということは、この新たなボクセルの各頂点のλ_２のみ残し新たなボクセルの各頂点以外の点のλ_２を削除することにより、解像度を下げることを意味している。このようにして作成された新たなボクセルに関し、上記と同様にしてλ_２＝−０．１の等値面（ポリゴン）を作成してそのポリゴンの数を集計する。ｘ，ｙ，ｚ方向に並ぶ、もともとのボクセルの数をそれぞれｍ，ｎ，ｏとしたとき、ボクセルの数をｍ／２，ｎ／２，ｏ／２；ｍ／３，ｎ／３，ｏ／３；…のように順次下げながら、すなわち、分解能を順次下げながら、ポリゴンの数が設定値以下となるまでポリゴンの作成を繰り返し、そのポリゴンの数の集計を行なう。

ここで、この設定値は、図１に示す可視化用計算機３０で渦の可視化のために十分な演算速度を持って処理することが可能なポリゴン数（例えば１００万個等）を意味している。この設定値は、可視化用計算機３０があらかじめ指定した値であってもよく、あるいは、一律に設定された、可視化用計算機として用い得る標準的なＰＣで十分な演算速度を持って処理できる値であってもよい。

設定値以下のポリゴン数が得られたら、そのポリゴン数が得られたボクセル１つずつについて、渦度の強さ｜ω｜が設定される。もともとの統合前の各ボクセルについて渦度の強さ｜ω｜が求められており（図６のステップＳ２６）、統合したボクセルについては、その統合したボクセルに含まれる複数の元々のボクセルの渦度の強さ｜ω｜の平均値や中央値等の代表値が設定される。

このようなポリゴン数削減処理（図６のステップＳ２９）により、流体空間内のλ_２や｜ω｜についての分解能は低下するが、可視化用計算機３０ではスムーズな表示が可能となる。このため渦の時間的変化（発生、消滅、分岐、合流等）の観察が容易となり、また、渦を別の角度から見るときの回転や、その他の渦の拡大、縮小、移動等の操作に対する応答性が向上する。また、表示の分解能は、渦抽出演算時点の分解能より低下するが、渦の中心線と領域とを表わすテキストデータの作成処理は分解能を下げる前の段階で行なっている。したがって、例えば分解能を下げることによってあらわれなくなるような、重要な過渡現象などもきちんと捉えることができる。

図６のデータ出力処理（Ｓ３０）では、図３のデータ出力部２７０に対応する処理が行なわれる。すなわちここでは、渦の中心を示す線分作成処理（Ｓ２８）で作成された渦の中心を示す線分を表わすテキストデータの出力処理が行なわれる。さらに、このデータ出力処理（Ｓ３０）では、ポリゴン数削減処理（Ｓ２９）で統合された後の新たなボクセルの各頂点に対応づけられたλ_２、および統合後の新たなボクセルに対応づけられた｜ω｜のデータも出力される。また、渦の中心線上のボクセル（値Ｖｃが対応づけられたボクセル）の位置座標も出力してもよい。ただし、本実施形態では、ポリゴン数削減処理（Ｓ２９）を実行する途中で作成されたポリゴンを表わすデータは出力されない。その理由は、ポリゴンを表わすデータはかなり大きなデータ量となり、通信時間や格納しておくためのメモリ容量上不利であること、および、必要であれば可視化用計算機３０で十分な演算速度でポリゴンを作成することが可能であること、である。

このようにして、図１に示す大規模計算機１０から出力されたデータは、ファイルサーバ２０に一旦格納され、可視化用計算機３０からの要求に応じて可視化用計算機３０に送られる。可視化用計算機３０では受け取ったデータに基づいて表示画面３１上にテキストデータに基づく渦の中心を示す線分やボクセルデータに基づく渦自体の表示が行なわれる。

ここでは、可視化用計算機３０における表示態様の如何を問うものではないが、以下では、可視化用計算機３０での典型的な表示態様を説明しておく。

最小限の表示としては、渦の中心を示す線分を表わすテキストデータに基づく渦の中心を示す線分の表示であり、渦の発生、消滅、合流、分岐、移動などを視認することができる。

詳細な表示としては、前述と同様にして、可視化用計算機３０側でＭａｒｃｈｉｎｇ Ｃｕｂｅ法を用いてポリゴンデータを作成し、表示画面３１上に渦の表面形状を表示することが挙げられる。

こうすると、テキストデータに基づく渦の中心を示す線分の表示では不明な、渦の詳細な形状を知ることができる。また、このようなポリゴンによる渦の表面形状の表示に加え、その表面のポリゴンを含むボクセルに対応づけられた｜ω｜をその大きさに応じた色で表示してもよい。例えば大きな値の｜ω｜については赤色、小さな値の｜ω｜については青色で表示する。こうすると、渦の強さが色分布で表示され、渦中心の方向も理解可能となる。

また、ポリゴンによる渦の表面形状の表示に加え、その表面形状を半透明で表示し、テキストデータに基づく渦の中心線や領域を重畳して表示してもよい。さらに、上記の各表示態様について、断面を表示して渦の内部を観察してもよい。

次に図面を参照しながら表示例を説明する。

図１７は、渦表示の一例を示す図である。

ここでは、立方体内に流体が満たされており、その立方体の上面を等速で移動させた場合を模擬して作成した流速分布を基にして抽出した渦を表わしたものである。その渦の表面の｜ω｜の大きさを色の濃淡であらわしてもよい。

図１８は、渦を、中心線が左右に延びる方向から見て半透明で表示した図である。また、図１９は、その渦を中心線が紙面に垂直な方向に延びる方向から見て、渦を半透明で示した図である。

図１８，図１９には中心線に相当する値Ｖｃが設定されたボクセルが示されている。

図２０，図２１は、それぞれ、図１８，図１９の表示画面上にさらに中心線と渦の中心を示す線分の領域を表わす線を重畳して示した図である。

本実施形態によれば、十分な描画速度をもって渦をこのような様々な形態で表示することが可能となる。また、本実施形態によれば、ボクセルを統合する前の高い分解能の演算で渦の中心を示す線分のテキストデータを作成しているため、重要な情報の欠落を避けることができる。

図２２は、本実施形態の変形例を示した図である。

ここには、図５に示す大規模計算機と同じ大規模計算機１０が示されている。ここで、図５では、流体空間内の流速分布を表わすデータはファイルサーバ２０（図１参照）にあらかじめ格納されていてそのデータを受け取るものとして説明した。これに対し、この図２２では、流体空間内の流速分布を求めるシミュレーションもこの大規模計算機１０で実行される。流体空間が領域分割され、各ＣＰＵ１１Ａ〜１１Ｄでは各領域を分担し、その分担した領域について先ずシミュレーションにより流速分布を算出し、次いで渦抽出を行なっている。渦抽出が終了すると、次の時刻についての流速分布算出のシミュレーションが行なわれ、その時刻の渦抽出が行なわれる。図２２の各ＣＰＵ１１Ａ〜１１Ｄでは、このようにして、流速分布算出のシミュレーションと渦抽出が交互に行なわれている。このように、本件では、渦抽出のみに特化している必要はなく、流速分布算出のシミュレーション用の計算機と兼ねてもよい。

１０ 大規模計算機
２０ ファイルサーバ
３０ 可視化用計算機
３１ 表示画面
４０ 通信回線網
１００ 渦抽出プログラム
１１０，２１０ データ取得部
１２０，２２０ 補間演算部
１３０，２３０ 第１の数値データ算出部
１４０，２４０ 第２の数値データ算出部
１５０，２５０ テキストデータ作成部
１６０，２６０ ポリゴン数削減処理部
１７０，２７０ データ出力部
２００ 渦抽出装置

Claims (10)

1. 流体空間内に分散した複数の各点の流速データを取得するデータ取得部と、
   前記データ取得部で取得した流速データに基づいて、前記流体空間内を複数のボクセルに分割したときの各ボクセルの各頂点ごとの、該各頂点が流体の渦領域の内側にあるか外側にあるかを指標する第１の数値データを算出する第１の数値データ算出部と、
   前記データ取得部で取得した流速データに基づいて流体の渦の中心線を求めるとともに、前記第１の数値データ算出部で算出された前記各頂点ごとの前記第１の数値データに基づいて該中心線の周りの渦領域の寸法を求め、該中心線と該寸法とを表わすテキストデータを作成するテキストデータ作成部と、
   前記第１の数値データ算出部で算出された前記第１の数値データの等値面で構成される、前記渦領域の表面を画定するポリゴンの数を計数し、前記複数のボクセルを、該ポリゴンの数が設定値以下となるまで前記流体空間が粗く分割された新たなボクセルに統合することにより、前記第１の数値データを、該新たなボクセルの各頂点に対応する第１の数値データにまで間引くポリゴン数削減処理部と、
   前記テキストデータ作成部で作成されたテキストデータおよび前記ポリゴン数削減処理部で得られた前記新たなボクセルの各頂点に対応する前記第１の数値データを出力するデータ出力部とを有することを特徴とする渦抽出装置。
2. 前記データ取得部で取得した流速データに基づいて、前記第１の数値データ算出部で算出された前記第１の数値データに基づいて判定される渦領域を構成する前記各ボクセルの渦の強さを指標する第２の数値データを算出する第２の数値データ算出部をさらに有し、
   前記ポリゴン数削減処理部は、前記新たなボクセルの各頂点に対応する第１の数値データにまで間引くとともに、前記第２の数値データ算出部で算出された第２の数値データに基づいて前記新たなボクセルそれぞれに対応する第２の数値データを作成するものであり、
   前記データ出力部は、前記テキストデータおよび前記第１の数値データに加え、さらに、前記新たなボクセルに対応する前記第２の数値データを出力するものであることを特徴とする請求項１記載の渦抽出装置。
3. 前記データ取得部で取得した流速データに基づく補間演算により、前記流体空間を構造格子に分割したときの複数のボクセルの各頂点ごとの流速データを算出する補間演算部を有することを特徴とする請求項１又は２記載の渦抽出装置。
4. 前記第１の数値データ算出部は、
   前記各ボクセルの各頂点の流速勾配行列を、代表的に
   

   

   とし、
   

   

   

   ただし、流体空間内のｘ，ｙ，ｚ方向の流速ｕ，ｖ，ｗをｑで代表させ、かつｘ，ｙ ，ｚをｒで代表させたとき、
   

   

   であり、またＪ^Ｔは行列Ｊの転置行列である。
   としたときの、行列Ａ＝Ｓ^２+Ｑ^２の３つの固有値λ_ａ，λ_ｂ，λ_ｃを値順にλ_１＜λ_２＜λ_３（但し、λ_１，λ_２，λ_３はそれぞれλ_ａ，λ_ｂ，λ_ｃのいずれか）と並べたときの中央の値の固有値λ_２を求め、該固有値λ_２に由来する数値データを前記第１の数値データとするものであることを特徴とする請求項１から３のうちいずれか１項記載の渦抽出装置。
5. 前記第２の数値データ算出部は、
   前記各ボクセルの流速ベクトルをＵとしたとき、
   

   

   但し、
   

   

   ×は外積
   ｜ω｜はωの大きさ
   を表わす。
   を算出し、該｜ω｜に由来する数値データを前記第２の数値データとするものであることを特徴とする請求項２記載の渦抽出装置。
6. 前記テキストデータ作成部は、流体空間内の関心点の周囲を取り巻く各点の流速データにより表わされる液体の流れの方向に基づいて該関心点が前記中心線上の点であるか否かを判定する処理を複数の関心点について実行することにより、該中心線を求めるものであることを特徴とする請求項１から５のうちいずれか１項記載の渦抽出装置。
7. 前記テキストデータ作成部は、前記渦領域の寸法として、前記中心線上の代表点から、該中心線が交わる方向に広がる面内での前記渦領域表面までの最短距離と最長距離とを求めるものであることを特徴とする請求項１から６のうちいずれか１項記載の渦抽出装置。
8. プログラムを実行する演算装置により実行され、該演算装置を、
   流体空間内に分散した複数の各点の流速データを取得するデータ取得部と、
   前記データ取得部で取得した流速データに基づいて、前記流体空間内を複数のボクセルに分割したときの各ボクセルの各頂点ごとの、該各頂点が流体の渦領域の内側にあるか外側にあるかを指標する第１の数値データを算出する第１の数値データ算出部と、
   前記データ取得部で取得した流速データに基づいて流体の渦の中心線を求めると共に、前記第１の数値データ算出部で算出された前記各頂点ごとの前記第１の数値データに基づいて該中心線の周りの渦領域の寸法を求め、該中心線と該寸法とを表わすテキストデータを作成するテキストデータ作成部と、
   前記第１の数値データ算出部で算出された前記第１の数値データの等値面で構成される、前記渦領域の表面を画定するポリゴンの数を計数し、前記複数のボクセルを、該ポリゴンの数が設定値以下となるまで前記流体空間が粗く分割された新たなボクセルに統合することにより、前記第１の数値データを、該新たなボクセルの各頂点に対応する第１の数値データにまで間引くポリゴン数削減処理部と、
   前記テキストデータ作成部で作成されたテキストデータおよび前記ポリゴン数削減処理部で得られた前記新たなボクセルの各頂点に対応する前記第１の数値データを出力するデータ出力部とを有する渦抽出装置として動作させることを特徴とする渦抽出プログラム。
9. 流体空間内に分散した複数の各点の流速データを取得するデータ取得過程と、
   前記データ取得過程で取得した流速データに基づいて、前記流体空間内を複数のボクセルに分割したときの各ボクセルの各頂点ごとの、該各頂点が流体の渦領域の内側にあるか外側にあるかを指標する第１の数値データを算出する第１の数値データ算出過程と、
   前記データ取得過程で取得した流速データに基づいて流体の渦の中心線を求めるとともに、前記第１の数値データ算出過程で算出された前記各頂点ごとの前記第１の数値データに基づいて該中心線の周りの渦領域の寸法を求め、該中心線と該寸法とを表わすテキストデータを作成するテキストデータ作成過程と、
   前記第１の数値データ算出過程で算出された前記第１の数値データの等値面で構成される、前記渦領域の表面を画定するポリゴンの数を計数し、前記複数のボクセルを、該ポリゴンの数が設定値以下となるまで前記流体空間が粗く分割された新たなボクセルに統合することにより、前記第１の数値データを、該新たなボクセルの各頂点に対応する第１の数値データにまで間引くポリゴン数削減処理過程と、
   前記テキストデータ作成過程で作成されたテキストデータおよび前記ポリゴン数削減処理過程で得られた前記新たなボクセルの各頂点に対応する前記第１の数値データを出力するデータ出力過程とをコンピュータが実行することを特徴とする渦抽出方法。
10. 流体空間内に分散した複数の各点の流速データに基づいて該流体空間内の渦を表わすデータを作成する渦抽出装置と、該渦抽出装置で作成されたデータを受け取って該渦を表示する渦表示装置とを有する渦抽出表示システムであって、
    前記渦抽出装置が、
    流体空間内に分散した複数の各点の流速データを取得するデータ取得部と、
    前記データ取得部で取得した流速データに基づいて、前記流体空間内を複数のボクセルに分割したときの各ボクセルの各頂点ごとの、該各頂点が流体の渦領域の内側にあるか外側にあるかを指標する第１の数値データを算出する第１の数値データ算出部と、
    前記データ取得部で取得した流速データに基づいて流体の渦の中心線を求めるとともに、前記第１の数値データ算出部で算出された前記各頂点ごとの前記第１の数値データに基づいて該中心線の周りの渦領域の寸法を求め、該中心線と該寸法とを表わすテキストデータを作成するテキストデータ作成部と、
    前記第１の数値データ算出部で算出された前記第１の数値データの等値面で構成される、前記渦領域の表面を画定するポリゴンの数を計数し、前記複数のボクセルを、該ポリゴンの数が設定値以下となるまで前記流体空間が粗く分割された新たなボクセルに統合することにより、前記第１の数値データを、該新たなボクセルの各頂点に対応する第１の数値データにまで間引くポリゴン数削減処理部と、
    前記テキストデータ作成部で作成されたテキストデータおよび前記ポリゴン数削減処理部で得られた前記新たなボクセルの各頂点に対応する前記第１の数値データを出力するデータ出力部とを有する渦抽出装置であることを特徴とする渦抽出表示システム。

Images