JP2018092473A - 流脈線可視化装置、流脈線可視化方法、および流脈線可視化プログラム - Google Patents

Abstract

【課題】構造体が変形しても流脈線の追跡を可能とする。【解決手段】流脈線可視化装置１０は、第１流脈線上の第１の位置５にある離散点を含む、解析空間１内の部分領域を、離散点についての解析対象領域４に決定する。次に流脈線可視化装置１０は、流体情報１１ｂに示される解析対象領域４内の流体の速度に基づいて、第２解析時刻における離散点上の粒子の移動先を示す第２の位置６を計算する。次に流脈線可視化装置１０は、構造体情報１１ａに示される解析対象領域４内に存在する構造体２の情報に基づいて、第２解析時刻における解析対象領域４内の構造体２の占有領域を特定する。次に流脈線可視化装置１０は、第１の位置５と第２の位置６とに基づいて、第２流脈線の占有領域への侵入の有無を判定する。そして流脈線可視化装置１０は、第２流脈線が占有領域に侵入しないとき、第２の位置６を通過する第２流脈線を表示する。【選択図】図１

Description

本発明は、流脈線可視化装置、流脈線可視化方法、および流脈線可視化プログラムに関する。

力学の一分野として流体力学がある。流体力学は、流体の振る舞いを記述する学問である。流体力学は、空気、水の流れはもとより、温度、濃度などの物理量の輸送が問題となる様々な産業分野に応用されている。例えば流体力学を用いることで、自動車の試作品の評価のための風洞実験に基づいて、自動車の空力特性を最適化することができる。ただし、風洞実験には多大なコストがかかる。そこで風洞実験に代えて、数値流体力学を用いることにより風洞実験を模擬したコンピュータシミュレーション（流体シミュレーション）が行われる。

コンピュータの性能向上に伴い、流体シミュレーションも飛躍的に発展している。その結果、流体シミュレーションは、航空機、自動車、鉄道車両、船舶などの空力特性の評価だけでなく、心臓、血管などの血流状態の解析にも応用されている。

流体シミュレーションを行った場合、解析結果を視覚により容易に把握できるように、解析結果が可視化される。流体シミュレーションの結果の可視化手段の一つとして、流脈線の表示がある。流脈線は、空間中のある点を通過した流体粒子を連ねてできる曲線である。風洞実験であれば、所定の場所から噴出された煙の軌跡が流脈線となる。すなわち、流体シミュレーションで流脈線を求めてその流脈線を表示することで、風洞実験における煙の軌跡と同様の流体中の粒子の動きを、風洞実験を行わずに可視化することができる。

流体シミュレーションに関する技術としては、例えば高速にシミュレーションを行い、流体の任意のシーンを短時間で滑らかにかつ細部の形状で表現する技術が考えられている。また、構造流体解析シミュレーションによる構造流体解析結果を、血管の異常診断に容易に適用可能とする技術もある。さらに、数値流体力学に精通していない医師などのユーザが適切な血流シミュレーションを行うことを可能とする装置も考えられている。その他、流体シミュレーションに関連する様々な論文が発表されている。

特開２００３−６５５２号公報 特開２０１５−９７７５９号公報 国際公開第２０１６／０５６６４２号

Tino Wienkauf and Holger Theisel, "Streak Lines as Tangent Curves of a Derived Vector Field", IEEE TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS, NOVEMBER/DECEMBER 2010, VOL. 16, NO. 6 Erwin Fehlberg, "NASA TECHNICAL REPORT", NASA TR R-315, JULY 1969 J.Donea, S. Huerta, J.-Ph. Ponthot and A. Rodrigues-Ferran, "Arbitrary Lagrangian-Eulerian methods", Encyclopedia of Computational Mechanics, John Wiley & Sons Ltd, November 2004, pp.413-437 Seiryo Sugiura, Takumi Washio, Asuka Hatano, Junichi Okada, Hiroshi Watanabe, Toshiaki Hisada, "Multi-scale simulations of cardiac electrophysiology and mechanics using the University of Tokyo heart simulator" ,Progress in Biophysics and Molecular Biology, October/November 2012, Volume 110, Pages 380-389

しかし、従来の流脈線の解析技術では、自動車などの周囲の空気の流れを解析する場合のように、解析空間中に存在する構造体が変形しないことが前提となっており、構造体が変形する場合は、流脈線を正確に追跡することが困難である。

１つの側面では、本件は、構造体が変形しても流脈線の追跡を可能とすることを目的とする。

１つの案では、流体シミュレーション上の時間内の複数の解析時刻における複数の粒子の連なりを示す流脈線を計算し、流脈線を表示する流脈線可視化装置が提供される。流脈線可視化装置は、記憶部と処理部とを有する。記憶部は、解析空間内の構造体の形状の時間変化を示す構造体情報と、解析空間内の流体が存在する領域内の複数の点における、流体の速度の空間変化又は／及び時間変化を示す流体情報とを記憶する。処理部は、以下の処理を実行する。

処理部は、第１解析時刻における第１流脈線に基づいて第２解析時刻における第２流脈線を計算する場合、第１流脈線上の第１の位置にある離散点を含む、解析空間内の部分領域を、離散点についての解析対象領域に決定する。次に処理部は、流体情報に示される解析対象領域内の流体の速度に基づいて、第２解析時刻における離散点上の粒子の移動先を示す第２の位置を計算する。次に処理部は、構造体情報に示される解析対象領域内に存在する構造体の情報に基づいて、第２解析時刻における解析対象領域内の構造体の占有領域を特定する。次に処理部は、第１の位置と第２の位置とに基づいて、第２流脈線の占有領域への侵入の有無を判定する。そして処理部は、第２流脈線が占有領域に侵入しないと判定したとき、第２の位置を通過する第２流脈線を表示する。

１態様によれば、構造体が変形しても流脈線の追跡が可能となる。

第１の実施の形態に係る流脈線可視化装置の構成の一例を示す図である。 第２の実施の形態のシステム構成例を示す図である。 可視化装置のハードウェアの一構成例を示す図である。 流脈線の計算例を示す図である。 可視化装置の機能を示すブロック図である。 弾性体情報ファイル群の一例を示す図である。 流体情報ファイル群の一例を示す図である。 流脈線可視化処理の手順の一例を示すフローチャートである。 時間発展計算処理の手順を示すフローチャートの前半である。 時間発展計算処理の手順を示すフローチャートの後半である。 流脈線のデータ例を示す図である。 グリッド点の位置が時間に対して連続的に変化する様子を示す図である。 時刻ｔ＝ｔ_kにおける場の情報を求める処理の手順の一例を示すフローチャートである。 時刻ｔ＝ｔ_kにおける弾性体場の情報を求める処理の手順の一例を示すフローチャートである。 計算の結果得られる可能性のある位置の例を示す図である。 状態変数の真理値表を示す図である。 状態判定の処理手順の一例を示すフローチャートである。 時間発展後の位置が流体内か否かの判定処理の手順の一例を示すフローチャートである。 移動ベクトルｄｒが通過する弾性体要素の検索処理の手順の一例を示すフローチャートである。 予測球半径の決定処理の手順の一例を示すフローチャートである。 時間分割による計算量削減の概念を示す図である。 時間分割数に応じた合計計算時間の変化の一例を示す図である。 要素の辺の長さの分布の一例を示す図である。 予測球半径に応じた計算コストの変化を示す図である。 移動距離の確率分布を示す図である。 確率分布を仮定した時の計算効率化曲線を示す図である。 流脈線の表示例を示す図である。 時間刻み幅を変えたときの軌跡の精度の変化を示す図である。

以下、本実施の形態について図面を参照して説明する。なお各実施の形態は、矛盾のない範囲で複数の実施の形態を組み合わせて実施することができる。
〔第１の実施の形態〕
まず、第１の実施の形態について説明する。第１の実施の形態は、構造体が変形しても流脈線の追跡が可能な流脈線表示装置である。

図１は、第１の実施の形態に係る流脈線可視化装置の構成の一例を示す図である。流脈線可視化装置１０は、流体シミュレーション上の時間内の複数の解析時刻における複数の粒子の連なりを示す流脈線を計算し、流脈線を表示する装置である。例えば流脈線可視化装置１０は、演算処理装置としてのプロセッサや主記憶装置としてのメモリを備えたコンピュータである。流脈線可視化装置１０は記憶部１１と処理部１２とを有する。

記憶部１１は、構造体情報１１ａと流体情報１１ｂとを記憶する。記憶部１１は、例えば流脈線可視化装置１０が有するメモリまたはストレージ装置である。
構造体情報１１ａは、解析空間１内の構造体２の形状の時間変化を示す情報である。流体情報１１ｂは、解析空間１内の流体が存在する領域（流体領域３ａ，３ｂ）内の複数の点における、流体の速度の時間変化を示す情報である。なお図１では流体領域３ａ，３ｂが個別に示されているが、図示していない部分で流体領域３ａと流体領域３ｂとは繋がっていてもいなくてもよいものとする。先天性心疾患などの症例では、左心と右心がつながっているものもある。構造体２は、例えば心臓である。この場合、流体は心臓内の血液である。構造体２が心臓であるとき、構造体情報１１ａと流体情報１１ｂとは、例えば心臓の拍動と血液の冠循環との連成解析を行うシミュレーションによって生成される。

処理部１２は、記憶部１１に格納されている構造体情報１１ａと流体情報１１ｂに基づいて、流脈線を計算し、計算した流脈線を表示する。処理部１２は、例えば流脈線可視化装置１０が有するプロセッサである。

処理部１２は、流脈線の可視化処理において、第１解析時刻（ｔ₀）における流脈線９ａ（第１流脈線）に基づいて第２解析時刻（ｔ₁）における流脈線９ｂ（第２流脈線）を計算するとき、以下のような処理を行う。

処理部１２は、第１解析時刻（ｔ₀）における流脈線９ａ上の第１の位置５にある離散点を含む、解析空間１内の部分領域を、離散点についての解析対象領域４に決定する（ステップＳ１）。解析対象領域４は、例えば第１の位置５を中心とする球状の領域である。この場合、処理部１２は、例えば、第１解析時刻（ｔ₀）と第２解析時刻（ｔ₁）とを含む時間帯における流体の最大速度を求め、第１解析時刻（ｔ₀）と第２解析時刻（ｔ₁）との差、および最大速度に基づいて、解析対象領域４の半径を算出する。なお、第１解析時刻（ｔ₀）と第２解析時刻（ｔ₁）との差は、流脈線計算における時間刻み幅である。

次に処理部１２は、流体情報に示される解析対象領域４内の流体の速度に基づいて、第２解析時刻（ｔ₁）における離散点上の粒子の移動先を示す第２の位置６を計算する（ステップＳ２）。具体的には、第２の位置６は、第１解析時刻（ｔ₀）において離散点に存在する仮想的な粒子が流体に流されたとき、第２解析時刻（ｔ₁）にどこに存在するのかを示している。

次に処理部１２は、構造体情報に示される解析対象領域４内に存在する構造体２の情報に基づいて、第２解析時刻（ｔ₁）における構造体２の解析対象領域４内の占有領域を特定する（ステップＳ３）。

次に処理部１２は、第１の位置５と第２の位置６とに基づいて、第２解析時刻（ｔ₁）における流脈線９ｂの占有領域への侵入の有無を判定する（ステップＳ４）。例えば処理部１２は、第２の位置６が流体内にあるか否かの第１の判定と、第１の位置５と第２の位置６とを結ぶ線分が構造体２の表面と交わるか否かの第２の判定とを行う。そして処理部１２は、第１の判定と第２の判定との結果に基づいて、流脈線９ｂの占有領域への侵入の有無を判定する。例えば処理部１２は、第２の位置６が流体内にあり、線分が構造体２の表面と交わっていなければ、流脈線９ｂは占有領域に侵入していないと判断する。

また処理部１２は、第２の位置６が流体内にあっても、線分が構造体２の表面と交わっている場合、流脈線９ｂは占有領域に侵入していると判断する。例えば図１の例において、第２の位置６が流体領域３ａにある場合、第２の位置６は流体内にあるものの、流脈線が、流体領域３ｂと流体領域３ａとの間にある構造体２の要素を横切ってしまう。従って、このような場合は、流脈線９ｂが構造体２の占有領域に侵入していると判断される。

さらに処理部１２は、第２の位置６が流体外にあるが、線分が構造体２の表面と交わってない場合、流脈線９ｂは占有領域に侵入していないと判断する。ただし、この場合、流脈線９ｂが流体外に飛び出してしまっているため、処理部１２は、流脈線９ｂを表示せずに計算を終了する。

処理部１２は、離散点が占有領域に侵入しないとき、第２の位置６を通過する流脈線９ｂを表示する（ステップＳ５）。図１の例では、表示画面７内に、構造体２の断面図上に重ねて、粒子発生源８を起点として、第１解析時刻（ｔ₀）における流脈線９ａと第２解析時刻（ｔ₁）における流脈線９ｂとが表示されている。

なお離散点が占有領域に侵入していると判定した場合、処理部１２は、第１解析時刻（ｔ₀）と第２解析時刻（ｔ₁）との間の時間帯内に１以上の第３解析時刻を設定し、第３解析時刻での離散点上の粒子の移動先を示す第３の位置を計算する。そして処理部１２は、第３の位置に基づいて第２の位置６を再計算する。再計算により離散点が占有領域に侵入していないと判断できた場合、処理部１２は、流脈線９ａ，９ｂを表示する。

このような流脈線の計算および表示を、所定時間内の複数の解析時刻それぞれについて実行することで、表示画面７には、形状が時間変化する構造体２の周囲の流体についての流脈線９ａ，９ｂを、正しく表示することができる。すなわち、時間変化する構造体２に流脈線が侵入してしまうような誤った流脈線を表示することが抑止されており、精度の高い流脈線が表示される。

しかも流脈線を計算する際に解析する構造体情報１１ａや流体情報１１ｂを、解析対象領域４内の情報に限定しているため、効率的に処理を実行することができる。このような効率的な流脈線の計算を可能としたため、構造体２の形状の変化を考慮した計算が可能となり、形状が時間変化する構造体２の周囲の流体についての流脈線の可視化が容易となっている。

なお処理部１２は、解析対象領域４の決定の際、例えば流体が存在する領域内の流体の速度が示された複数の点の間隔に基づいて、解析対象領域４の半径の最小値を設定してもよい。この場合、処理部１２は、算出した半径が最小値より小さい場合、最小値を解析対象領域４の半径に決定する。例えば、複数の点の内の隣接する２点間の間隔の最大値を、解析対象領域４の半径の最小値とすることができる。これにより、解析対象領域４の半径を小さくしすぎて、第２の位置６を計算できるだけの、解析対象領域４内の流体の速度の情報が得られないような事態を抑止することができる。

なお処理部１２は、第１解析時刻（ｔ₀）と第２解析時刻（ｔ₁）との差を、第１解析時刻（ｔ₀）と第２解析時刻（ｔ₁）とを含む時間帯における流体の最大速度に乗算した値よりも、解析対象領域４の半径が小さいことを許容してもよい。この乗算の結果は、流体内の最も速度が早い位置にある離散点の移動距離（最大移動距離）である。最大移動距離を、解析対象領域４の半径の最大値とすることができる。例えば、処理部１２は、解析対象領域４の半径の最小値を、最大移動距離よりも小さい値に設定する。この場合、離散点上の粒子の移動先が解析対象領域４外になることがあり得る。そこで処理部１２は、第２の位置６の計算において、離散点上の粒子の移動先が解析対象領域４外となった場合、解析対象領域４を拡大して第２の位置６を再計算する。例えば処理部１２は、解析対象領域４を、最大移動距離を与える半径の球領域として再計算する。

このように解析対象領域４の半径が、離散点の最大移動距離よりも小さいことを許容することで、解析対象領域４の半径を小さくすることができる。解析対象領域４の半径が小さくなれば、解析対象の情報量も少なくなり、処理が効率化される。ただし、決定した半径の解析対象領域４内の第２の位置６の計算が失敗した場合、代わりに解析対象領域４の半径を最大値（最大移動距離）にして再計算することとなり、余分な再処理が発生する。そこで、処理部１２は、失敗するまでは、失敗した場合の第２の位置６の計算の量も含めて計算量を最小化する半径を算出し、算出した半径が最小値より小さい場合、最小値を解析対象領域４の半径に決定することができる。この場合、処理部１２は、例えば、離散点上の粒子の移動先が解析対象領域４外となる確率に基づいて、最も処理効率がよい半径を決定する。すなわち処理部１２は、離散点上の粒子の移動先が解析対象領域４外となる場合の処理の増加量と、解析対象領域４をより小さくすることによる処理の減少量とに基づいて、最も処理量が少なくなるように、解析対象領域４の半径を決定する。これにより、処理効率が向上する。

しかも、解析対象領域４の半径を小さくすることが可能となることで、有限要素モデルの良し悪しによって、一部の品質の悪いメッシュによって解析対象領域４の半径の最小値が制限されることを抑止できる。その結果、品質の悪いメッシュが存在することによる計算量の増大を抑止できる。

なお上記の説明では、流脈線９ｂの占有領域への侵入の有無を判定例として、第１の判定と第２の判定との両方を用いる例を示しているが、簡易的には、第２の判定のみを用いることもできる。すなわち、第１の位置５と第２の位置６とを結ぶ線分が構造体２の表面と交わるか否かの第２の判定のみで、流脈線９ｂの構造体２の占有領域への侵入の有無を判定することもできる。この場合、処理部１２は、線分が構造体２の表面と少なくとも１回交わっていれば、流脈線９ｂが構造体２の占有領域へ侵入すると判断する。例えば流体領域３ａ，３ｂが閉空間であり、心臓の動脈のような流体の流出先がない場合であれば、第２の判定のみでも正確な判定となる。

〔第２の実施の形態〕
次に第２の実施の形態について説明する。第２の実施の形態は、心臓の動きと共に、心臓内の血流の流脈線を可視化する可視化装置である。

例えば数値流体解析を用いれば、心臓中の血流輸送など技術的、倫理的観点から測定が困難である系についても流体の振る舞いを模擬することが可能である。そのため、心臓中での血流輸送の不全が問題となる先天性心疾患などでは治療法の検討に数値流体解析が使用されており、数値流体解析は重要な技術である。そして可視化装置により、このような数値流体解析の結果を可視化することで、医者などの医療関係者が解析結果を容易に把握することができ、治療方針を立てやすくなる。

図２は、第２の実施の形態のシステム構成例を示す図である。可視化装置１００は、ネットワーク２０を介して心臓シミュレータ２００に接続されている。心臓シミュレータ２００は、心筋の運動と冠循環とのシミュレーションを行うコンピュータである。可視化装置１００は、心臓シミュレータ２００からシミュレーション結果を取得する。そして可視化装置１００は、シミュレーション結果に基づいて流脈線を計算し、計算した流脈線を表示する。シミュレーション結果には、例えば各時刻における心臓の形状を示す３次元モデル、血管中の血液の速度、心筋または血液の物性値などの情報が含まれる。

図３は、可視化装置のハードウェアの一構成例を示す図である。可視化装置１００は、プロセッサ１０１によって装置全体が制御されている。プロセッサ１０１には、バス１０９を介してメモリ１０２と複数の周辺機器が接続されている。プロセッサ１０１は、マルチプロセッサであってもよい。プロセッサ１０１は、例えばＣＰＵ（Central Processing Unit）、ＭＰＵ（Micro Processing Unit）、またはＤＳＰ（Digital Signal Processor）などの演算処理装置である。プロセッサ１０１がプログラムを実行することで実現する機能の少なくとも一部を、ＡＳＩＣ（Application Specific Integrated Circuit）、ＰＬＤ（Programmable Logic Device）などの電子回路で実現してもよい。

メモリ１０２は、可視化装置１００の主記憶装置として使用される。メモリ１０２には、プロセッサ１０１に実行させるＯＳ（Operating System）のプログラムやアプリケーションプログラムの少なくとも一部が一時的に格納される。また、メモリ１０２には、プロセッサ１０１による処理に必要な各種データが格納される。メモリ１０２としては、例えばＲＡＭ（Random Access Memory）などの揮発性の半導体記憶装置が使用される。

バス１０９に接続されている周辺機器としては、ストレージ装置１０３、グラフィック処理装置１０４、入力インタフェース１０５、光学ドライブ装置１０６、機器接続インタフェース１０７およびネットワークインタフェース１０８がある。

ストレージ装置１０３は、内蔵した記録媒体に対して、電気的または磁気的にデータの書き込みおよび読み出しを行う。ストレージ装置１０３は、コンピュータの補助記憶装置として使用される。ストレージ装置１０３には、ＯＳのプログラム、アプリケーションプログラム、および各種データが格納される。なお、ストレージ装置１０３としては、例えばＨＤＤ（Hard Disk Drive）やＳＳＤ（Solid State Drive）を使用することができる。

グラフィック処理装置１０４には、モニタ２１が接続されている。グラフィック処理装置１０４は、プロセッサ１０１からの命令に従って、画像をモニタ２１の画面に表示させる。モニタ２１としては、ＣＲＴ（Cathode Ray Tube）を用いた表示装置や液晶表示装置などがある。

入力インタフェース１０５には、キーボード２２とマウス２３とが接続されている。入力インタフェース１０５は、キーボード２２やマウス２３から送られてくる信号をプロセッサ１０１に送信する。なお、マウス２３は、ポインティングデバイスの一例であり、他のポインティングデバイスを使用することもできる。他のポインティングデバイスとしては、タッチパネル、タブレット、タッチパッド、トラックボールなどがある。

光学ドライブ装置１０６は、レーザ光などを利用して、光ディスク２４に記録されたデータの読み取りを行う。光ディスク２４は、光の反射によって読み取り可能なようにデータが記録された可搬型の記録媒体である。光ディスク２４には、ＤＶＤ（Digital Versatile Disc）、ＤＶＤ−ＲＡＭ、ＣＤ−ＲＯＭ（Compact Disc Read Only Memory）、ＣＤ−Ｒ（Recordable）／ＲＷ（ReWritable）などがある。

機器接続インタフェース１０７は、可視化装置１００に周辺機器を接続するための通信インタフェースである。例えば機器接続インタフェース１０７には、メモリ装置２５やメモリリーダライタ２６を接続することができる。メモリ装置２５は、機器接続インタフェース１０７との通信機能を搭載した記録媒体である。メモリリーダライタ２６は、メモリカード２７へのデータの書き込み、またはメモリカード２７からのデータの読み出しを行う装置である。メモリカード２７は、カード型の記録媒体である。

ネットワークインタフェース１０８は、ネットワーク２０に接続されている。ネットワークインタフェース１０８は、ネットワーク２０を介して、他のコンピュータまたは通信機器との間でデータの送受信を行う。

以上のようなハードウェア構成によって、第２の実施の形態の処理機能を実現することができる。なお、第１の実施の形態に示した装置も、図３に示した可視化装置１００と同様のハードウェアにより実現することができる。

可視化装置１００は、例えばコンピュータ読み取り可能な記録媒体に記録されたプログラムを実行することにより、第２の実施の形態の処理機能を実現する。可視化装置１００に実行させる処理内容を記述したプログラムは、様々な記録媒体に記録しておくことができる。例えば、可視化装置１００に実行させるプログラムをストレージ装置１０３に格納しておくことができる。プロセッサ１０１は、ストレージ装置１０３内のプログラムの少なくとも一部をメモリ１０２にロードし、プログラムを実行する。また可視化装置１００に実行させるプログラムを、光ディスク２４、メモリ装置２５、メモリカード２７などの可搬型記録媒体に記録しておくこともできる。可搬型記録媒体に格納されたプログラムは、例えばプロセッサ１０１からの制御により、ストレージ装置１０３にインストールされた後、実行可能となる。またプロセッサ１０１が、可搬型記録媒体から直接プログラムを読み出して実行することもできる。

次に流脈線について説明する。
図４は、流脈線の計算例を示す図である。可視化装置１００は、解析対象の空間中に粒子発生源３０を定義する。解析が始まると、粒子群３３が粒子発生源３０から連続的に射出される。粒子群３３は流れ場が時間に対して不変であれば、一定の曲線（流線）になる。流れ場が時間依存して変化すると、粒子群３３が作り出す曲線も時々刻々と変化していく。この曲線が流脈線３１，３２である。図４では、時刻ｔ₀における粒子群３３の連なりを示す流脈線３１と時刻ｔ₁における粒子群３３の連なりを示す流脈線３２とが示されている。流脈線３１，３２は、障害物３５の影響で複雑な曲線となる。

流脈線３１，３２は時間変動する流れ場に沿って、どのように粒子群３３が輸送されていくかを可視化する有力な方法である。例えば障害物３５を自動車とした場合を考える。開発した自動車の空気抵抗を可視化するため、自動車前方から空気を送ると共に、自動車の前方に粒子発生源３０を配置する。そして可視化装置１００内での流体シミュレーションにおいて、粒子群３３を連続的に射出し、粒子群３３の軌跡が流脈線３１，３２として測定される。流脈線３１，３２は流体の輸送を直接記述し、可視化できるため、様々な分野に応用可能である。

流脈線の計算や、可視化の研究は多数行われており、乱流や不安定流などを対象とした研究は多数ある。しかし、心臓のように流体から見て壁面になるような弾性体が大きく変形を起こす場合の流脈線の可視化については先行研究があまりない。心臓は周期的に拍動し、拡張と収縮を繰り返すため、大変形を起こす系の典型例である。そして、心臓のポンプ作用はこの周期的な運動が重要な役割を果たしているため、弾性体が周期的に大変形を起こす系における血流の輸送を評価することは、心疾患の治療法を考慮するうえで重要になる。

生体シミュレーションの研究では、心臓の挙動をコンピュータ上でシミュレートすることも行われている。コンピュータ上でシミュレーションを行うことで、実際に手術を実行しなくても、手術を実行した場合の治療効果を評価できる。そのため生体シミュレーションを用いれば、医師は実際に手術をする前に、最もよい治療法を検討することが可能である。特に、心臓シミュレーションは、３次元的に複雑な構造を持つ心臓が対象である。そして、心臓の挙動も動的に変化する。心臓の挙動に合わせて、心臓内の血流の輸送を表す流脈線が可視化できれば、心臓の状態を視覚的に容易に把握することができる。視覚的に分かりやすく表示することは、判断の間違いを防ぐためにも有用である。

心臓内の血流の流脈線を可視化する際に障害となるのは以下の点である。
１．心筋（弾性体）が大きく変形するとき、心筋付近で流跡線、流脈線の挙動を正確に追跡することが困難なことである。

２．流跡線が常に１点のみ計算すればよいのに対して、流脈線は線分を構成するＮ点全てを計算しなければならず、計算量が非常に多いことである。
３．有限要素モデルの一部の品質の悪いメッシュが、全体の計算量を増加させてしまうことである。

そこで、第２の実施の形態の可視化装置１００は、以下の機能により、実現可能な計算量で高精度の流脈線を可視化できるようにしている。
１−１：可視化装置１００は、流脈線の各点が可動領域外の心筋に入り込んだり、シミュレーション対象の系の外に出たりしたかどうかを正確に判定する。

１−２：可視化装置１００は、流脈線上の点が可動領域外にある場合、流線の微分方程式のパラメータの時間刻み幅を調節し、可動領域外に出ないように制御する。
１−３：可視化装置１００は、任意の時刻における場の情報を推定するため、補間法を用いて場を補間する。

２−１：１−３の機能を適用すると計算量が増加するため、可視化装置１００は、流脈線の点が動く最大距離を評価し、その距離を半径とする予測球内部の場の情報のみを評価することでデータ容量に関係なく計算量を一定にする。

２−２：可視化装置１００は、時間刻み幅を分割することで予測球半径を小さくし、予測球を使用しない場合より計算量を減らすと共に、同時に精度も向上させる。
３−１：大部分の計算は品質のよいメッシュに対して行われるのを利用し、可視化装置１００は、品質のよいメッシュと決めて投機的に計算を実行し、失敗した場合には高精度計算を行うことで、計算量を削減する。投機的な計算とは、流脈線を構成する点の移動先が予測球の外部になる可能性があることを許容して、予測球の半径を小さくすることである。計算に失敗する場合とは、流脈線を構成する点の移動先が予測球内に存在しない場合である。

３−２：３−１にて投機的に実行するため、可視化装置１００は確率モデルを用意し、失敗時のペナルティまで含めて計算量が最小になるようなパラメータセットを決定する。
これらの機能が可視化装置１００に実装されていることで、以下の効果が得られている。

１．心筋（弾性体）が大きく変形をする場合にも心筋（弾性体）の運動を考慮して流脈線を計算できる。
２．流脈線の各点を、予測球を用いて高速・高精度に計算ができる。

３．確率モデルを用い、計算コストを最小化する予測球半径の決定ができる。
以下、可視化装置１００の機能を詳細に説明する。
図５は、可視化装置の機能を示すブロック図である。可視化装置１００は、記憶部１１０を有する。記憶部１１０は、心臓シミュレータ２００から取得したシミュレーション結果を記憶する。例えば数値流体力学シミュレーションが心臓シミュレータ２００で実施され、動的に変化する弾性体と流体の場のシミュレーション結果がＬ個（Ｌは１以上の整数）の時刻ｔ₀，ｔ₁，・・・，ｔ_Lでファイルに保存される。例えば心筋の情報と血流の情報が別々のファイルとして、記憶部１１０に保存される。図５の例では、時刻ごとの心筋の情報が弾性体情報ファイル群１１１として保存され、時刻ごとの血流の情報が流体情報ファイル群１１２として保存されている。

このシミュレーション結果を分析することで、血流の輸送情報を記述する流脈線を計算することができる。なお、シミュレーション結果として出力される時間間隔Δｔ_i＝ｔ_i+1−ｔ_iは心臓シミュレータ２００が微分方程式を解く場合の時間間隔と一致していなくてもよい。情報量を削減するため、間引かれたシミュレーション結果が出力されることが一般的である。そこで可視化装置１００は、流脈線を高精度に求めるため、複数の時刻での出力ファイルを用いてターゲットとなる時刻での種々の物理量を、補間法などを用いて推定する。

次に可視化装置１００の処理機能について説明する。可視化装置１００は、情報読み込み部１２０、流脈線計算部１３０、表示処理部１４０、および解析部１５０を有する。
情報読み込み部１２０は、流体解析の結果を示すファイルを記憶部１１０から読み込む。流脈線計算部１３０は、読み込んだ情報を用いて流脈線を計算する。表示処理部１４０は、得られた結果を可視化する。

解析部１５０は、情報読み込み部１２０、流脈線計算部１３０、および表示処理部１４０が共通に使用する機能をまとめたものである。情報読み込み部１２０、流脈線計算部１３０、および表示処理部１４０は、具体的な解析処理を実行する際に、解析部１５０に実行する処理を依頼して、結果を得る。

解析部１５０は、流体情報解析部１５１、弾性体情報解析部１５２、予測球半径計算部１５３、予測球内部計算部１５４、移動境界衝突判定部１５５、および表示様式決定部１５６を有する。流体情報解析部１５１は、流体の速度場や離散点の位置、境界面の解析などを行う。弾性体情報解析部１５２は、心筋など流体ではない弾性体の離散点の位置や境界面の解析を行う。予測球半径計算部１５３は、流脈線の計算時の計算速度と計算精度の向上のために用いる予測球の半径を決定する。予測球内部計算部１５４は、予測球内部の速度場、心筋位置などを計算する。移動境界衝突判定部１５５は、流脈線上の点が計算上の誤差のため、心筋内部に入り込んだかどうかを判定する。表示様式決定部１５６は、得られた流脈線の表示の仕方を決定する。

なお、図５に示した各要素の機能は、例えば、その要素に対応するプログラムモジュールをコンピュータに実行させることで実現することができる。
次に、シミュレーション結果として得られる情報について詳細に説明する。

図６は、弾性体情報ファイル群の一例を示す図である。弾性体情報ファイル群１１１は、シミュレーションの時刻ごとの弾性体情報ファイル１１１ａ，１１１ｂ，・・・の集合である。各弾性体情報ファイル１１１ａ，１１１ｂ，・・・には、例えば「stru(X).inp」というファイル名が付与される。ファイル名中の「Ｘ」は、シミュレーションの時刻が古い順に付与された昇順の番号である。

弾性体情報ファイル１１１ａ，１１１ｂ，・・・内には、対応する時刻における心臓の形状を表す心筋データが含まれている。心筋データには、グリッド（３次元空間内に配置された頂点）のｘ、ｙ、ｚの座標値、心筋に含まれる４面体の要素（ＴＥＴＲＡ）の４つの頂点を示すグリッドのＩＤ、各要素に加わる力などの情報が含まれる。

図７は、流体情報ファイル群の一例を示す図である。流体情報ファイル群１１２は、シミュレーションの時刻ごとの流体情報ファイル１１２ａ，１１２ｂ，・・・の集合である。各流体情報ファイル１１２ａ，１１２ｂ，・・・には、例えば「flui(Y).inp」というファイル名が付与される。ファイル名中の「Ｙ」は、シミュレーションの時刻が古い順に付与された昇順の番号である。

流体情報ファイル１１２ａ，１１２ｂ，・・・内には、対応する時刻における血液の流れを表す血流データが含まれている。血流データには、グリッド（３次元空間内に配置された頂点）のｘ、ｙ、ｚの座標値、血管内の４面体の要素（ＴＥＴＲＡ）の４つの頂点を示すグリッドのＩＤ、各要素内の血液の流れる方向と早さを示す速度場ベクトルが含まれる。

可視化装置１００は、図６と図７に示したシミュレーション結果に基づいて、流脈線を計算し可視化する。以下、流脈線の可視化処理について詳細に説明する。
図８は、流脈線可視化処理の手順の一例を示すフローチャートである。以下、図８に示す処理をステップ番号に沿って説明する。

まず、流脈線計算部１３０は、ステップＳ１０１〜Ｓ１０４において、流脈線の初期設定を行う。
［ステップＳ１０１］流脈線計算部１３０は、計算する流脈線の本数Ｍ（Ｍは１以上の整数）を決定する。例えば流脈線計算部１３０は、ユーザから入力された値を、流脈線の本数Ｍとする。

［ステップＳ１０２］流脈線計算部１３０は、流脈線の計算回数Ｎ（Ｎは１以上の整数）を決定する。以下、流脈線の計算回数Ｎを「時間発展回数」と呼ぶ。例えば流脈線計算部１３０は、ユーザから入力された値を、時間発展回数Ｎとする。

流脈線は時間と共に変化するので、時間発展回数Ｎを決定することで、流脈線の結果が出力される時刻列ｔ₀，ｔ₁，・・・，ｔ_Nが決定される。指定された時間発展回数Ｎがファイルの個数Ｌ（Ｌは１以上の整数）より大きい場合は、流脈線計算部１３０は、Ｌ＋１番目を心臓の２周期目の拍動として扱い、ｔ₀でのファイルを使用する。時刻列は、例えば０．０１秒間隔の等間隔に設定される。ただし、時刻列の時間間隔は、等間隔でなくてもよい。

［ステップＳ１０３］流脈線計算部１３０は、流脈線の発生源の座標を決定する。例えば流脈線計算部１３０は、ユーザから指定された解析空間上の点を、発生源の座標とする。ユーザは、例えば心筋情報と血流情報を参照しながら、空間の一点を指定する。流脈線計算部１３０は、指定された点の座標を、座標ベクトルＸ₀として読み込む。流脈線の本数が１本の場合は、流脈線の発生源は座標ベクトルＸ₀に設定される。流脈線の本数が複数本ある場合、流脈線計算部１３０は、座標ベクトルＸ₀を中心とした半径ｒ（ｒは正の実数）の球の中から流脈線ｌ_jの発生源をランダムに決定する。なお発生源は、血流中の座標から選ばれる。流脈線計算部１３０は、決定した発生源の座標ベクトルＸ_jを流脈線発生源として設定する。

その後、流脈線計算部１３０は、ｊ番目（ｊ＝１，２，・・・，Ｍ）の流脈線ｌ_jの初期設定を以下のようにして行う。
ｊ番目の流脈線ｌ_jは時間発展回数Ｎと同じ分だけの離散点で構成される。そこで流脈線計算部１３０は、流脈線ｌ_jに含まれる各離散を示す点Ｐ_ij（ｉ＝０，１，２，・・・，Ｎ）を生成する。流脈線計算部１３０は、各離散点の初期値の座標を、粒子発生源の座標とする。

なお、時刻ｔ＝ｔ_iにおける流脈線の計算では、点Ｐ_ij（ｉ＝０，１，２，・・・，ｉ）が、発生源から射出された流脈線粒子の位置として、時間発展計算の対象になる。点Ｐ_ij（ｉ＝ｉ＋１，・・・，Ｎ）に対応する流脈線粒子はまだ発生源から射出されていないと考えるため、時刻ｔ＝ｔ_iにおける流脈線の計算では計算対象とはならない。そして流脈線計算部１３０は、離散点のうちｉの値が小さい順に計算をする。そのため、先に計算される離散点ほど、粒子発生源から射出されてからの時間は長くなる。

［ステップＳ１０４］流脈線計算部１３０は、流脈線の点が大動脈など、系の流体境界から出た場合のための設定を行う。流脈線上の点Ｐ_ijが大動脈など流体境界から系外に出た場合、系外では流体の速度場が定義されていないため、次の時刻でのＰ_ijの計算を実行することはできない。そこで流脈線計算部１３０は、すべての点Ｐ_ijに領域内判定フラグを、各離散点のパラメータとして設定する。領域内判定フラグは、領域内にある場合を「Ｔ」とし、大動脈など流体の流れに従い領域外に出た場合を「Ｆ」とする。初期状態では流体内部に点Ｐ_ijが存在しているため、流脈線計算部１３０は、すべての離散点の領域判定フラグにＴをセットする。

［ステップＳ１０５］流脈線計算部１３０は、インデックスｉ＝１から順番に、インデックスｉ＝１，２，・・・，Ｎ−１について、ステップＳ１０６〜Ｓ１１０の処理を繰り返し実行する。

［ステップＳ１０６］流脈線計算部１３０は、インデックスｊ＝１から順番に、インデックスｊ＝１，２，・・・，Ｍについて、ステップＳ１０７〜Ｓ１０９の処理を繰り返し実行する。

［ステップＳ１０７］流脈線計算部１３０は、時間発展の起点になる時刻を決定し、メモリ１０２に保持する。計算開始時刻は第ｉ回目の計算ではｔ＝ｔ_iに設定される。時間発展の終端時刻はｔ＝ｔ_i+1に設定される。

［ステップＳ１０８］流脈線計算部１３０は、時間ｔ_i≦ｔ≦ｔ_i+1における時間発展計算を行う。時間発展計算により、各時刻ｔ＝ｔ_iにおいて、流脈線発生源から射出されたすべての点Ｐ_ij（ｉ＝０，１，２，・・・，ｉ）が時間発展の対象となり、線分上のすべての点が時々刻々と更新されていく。流脈線ｌ_j上の点Ｐ_ijの時間［ｔ_i，ｔ_i+1］の時間発展計算の結果、座標値が返される。これが次の時刻ｔ＝ｔ_i+1での座標Ｐ_i+1,jとして設定される。

［ステップＳ１０９］流脈線計算部１３０は、得られた計算結果をメモリ上に保存する。計算結果は、表示処理部１４０によって流脈線が描画される。また流脈線計算部１３０は、ファイルに流脈線の座標値を出力することもできる。

［ステップＳ１１０］流脈線計算部１３０は、ステップＳ１０７〜Ｓ１０９の処理を１回実行するごとにインデックスｊに１を加算し、インデックスｊ＝Ｍの処理が終了したら、処理をステップＳ１１１に進める。

［ステップＳ１１１］流脈線計算部１３０は、ステップＳ１０６〜Ｓ１１０の処理を１回実行するごとにインデックスｉに１を加算し、インデックスｊ＝Ｎ−１の処理が終了したら、流脈線可視化処理を終了する。

次に、時間発展計算処理（ステップＳ１０８）の詳細について説明する。
図９は、時間発展計算処理の手順を示すフローチャートの前半である。以下、図９に示す処理をステップ番号に沿って説明する。

［ステップＳ１２１］流脈線計算部１３０は、点Ｐ_ijの領域内判定フラグを読み込む。そして流脈線計算部１３０は、領域内判定フラグが「Ｔ」か否かを判断する。領域内判定フラグが「Ｔ」であれば、流脈線計算部１３０は、処理をステップＳ１２２に進める。領域内判定フラグが「Ｆ」であれば、流脈線計算部１３０は、この点は既に領域外に出たものとして扱い、計算を行わず、時間発展計算処理を終了する。これにより、領域内判定フラグが「Ｆ」の場合、座標値は更新されないこととなる。

［ステップＳ１２２］流脈線計算部１３０は、点Ｐ_ijの座標値を読み込む。
［ステップＳ１２３］流脈線計算部１３０は、計算開始時刻ｔ＝ｔ_iおよび計算終了時刻ｔ＝ｔ_i+1を、メモリ１０２から読み込む。

［ステップＳ１２４］流脈線計算部１３０は、計算開始時刻ｔ＝ｔ_iに対応する、流体部分のグリッド情報ベクトルＢ（ベクトルｒ，ｔ_i）と速度場ベクトルｖ（ベクトルｒ，ｔ_i）を、情報読み込み部１２０を介してファイルflui(i).inpから読み込む。また流脈線計算部１３０は、心筋部分の弾性体のグリッド情報（心筋の構造情報）であるベクトルＭ（ベクトルｒ，ｔ_i）を、情報読み込み部１２０を介してファイルstru(i).inpから読み込む。

［ステップＳ１２５］流脈線計算部１３０は、計算終了時刻ｔ＝ｔ_i+1における流体部分のグリッド情報ベクトルＢ（ベクトルｒ，ｔ_i+1）と速度場ベクトルｖ（ベクトルｒ，ｔ_i+1）を、情報読み込み部１２０を介してファイルflui(i+1).inpから読み込む。また流脈線計算部１３０は、ベクトルＭ（ベクトルｒ，ｔ_i+1）を、情報読み込み部１２０を介してファイルstru(i+1).inpから読み込む。

［ステップＳ１２６］流脈線計算部１３０は、flui(i).inpとflui(i+1).inpの各グリッド点の速度場ベクトルのノルム（速度場ベクトルの長さ）を計算し、補間式を用いることにより当該時間区間における速度の最大値を求め、求めた最大値をＶ_maxとしてメモリ１０２に保存する。

［ステップＳ１２７］区間［ｔ_i，ｔ_i+1］をそのまま１回の時間発展で計算すると精度が不十分である。そのため流脈線計算部１３０は、区間［ｔ_i，ｔ_i+1］をＮ_div（Ｎ_divは１以上の整数）等分し、中間時刻（ｔ＝ｔ_k（ｋ＝１，２，．．．，Ｎ_div））を決定する。これにより、区間を［ｔ_i，ｔ_i＋Δｔ］，［ｔ_i＋Δｔ，ｔ_i＋２Δｔ］，・・・，［ｔ_i＋（Ｎ_div−１）Δｔ，ｔ_i+1］に分割して解くこととなる。分割数Ｎ_divは、流脈線計算部１３０が最適値を自動で決定する。なお、分割数Ｎ_divを外部から与えることもできる。その後、処理がステップＳ１３１（図１０参照）に進められる。

図１０は、時間発展計算処理の手順を示すフローチャートの後半である。以下、図１０に示す処理をステップ番号に沿って説明する。
［ステップＳ１３１］流脈線計算部１３０は、中間時刻（ｔ＝ｔ_k（ｋ＝１，２，．．．，Ｎ_div））の、ｋ＝１からｋ＝Ｎ_divまで、ステップＳ１３２〜Ｓ１４０の処理を繰り返し実行することで、時間発展計算を行う。これにより、各中間時刻ｔ＝ｔ_kにおける座標ベクトルｒ_kが得られ、時刻ｔ＝ｔ_i+1における点Ｐ_i+1jの座標ベクトルｒ_Ndiv＝ベクトルｒ_i+1が得られる。

［ステップＳ１３２］流脈線計算部１３０は、時刻ｔ＝ｔ_kにおける場の情報（流体の構造情報ベクトルＢ（ベクトルｒ，ｔ_k）と速度場ベクトルｖ（ベクトルｒ，ｔ_k））を求める。この処理の詳細は後述する（図１３参照）。

［ステップＳ１３３］流脈線計算部１３０は、時刻ｔ＝ｔ_kにおける弾性体場の情報（心筋の構造情報ベクトルＢ（ベクトルｒ，ｔ_k））を求める。この処理の詳細は後述する（図１４参照）。

［ステップＳ１３４］流脈線計算部１３０は、時間ｄｔ＝ｔ_k+1−ｔ_k分だけ時間発展を行い、ベクトルｒ_k+1を求める。
ステップＳ１３２〜Ｓ１３４における各中間時刻での時間発展は以下のようにして行うことができる。点Ｐ_i+1jが時刻ｔ＝ｔ_kで、ベクトルｒ_kにあったとすると、流脈線の方程式は式（１）で与えられる。

ここでベクトルｖ（ベクトルｒ，ｔ）は時刻ｔ、位置ベクトルｒでの速度（場）である。従って、時間Δｔ後の座標は常微分方程式である式（１）を数値的に解くことによって得られる。式（１）を４次ルンゲ・クッタ法で解くと以下のような計算式が得られる。

時刻ｔ＝ｔ_kでの任意の位置ベクトルｒでの速度場ベクトルｖ（ベクトルｒ，ｔ_k）は、ベクトルｖ（ベクトルｒ，ｔ_i）、ベクトルｖ（ベクトルｒ，ｔ_i+1）、ベクトルＢ（ベクトルｒ，ｔ_i）、ベクトルＢ（ベクトルｒ，ｔ_i+1）、ベクトルＭ（ベクトルｒ，ｔ_i）、および時刻ｔ＝ｔ_kでの移動境界面Ｓ_kの情報から計算することができる。以下、移動境界面を単に「境界面」と呼ぶ。時刻ｔ＝ｔ_kでの境界面Ｓ_kの情報は、ベクトルＭ（ベクトルｒ，ｔ_i）から求められる。従って、式（３）、式（４）、式（５）、式（６）を計算することができる。これらの結果を式（２）に代入することで、ベクトルｒ_k+1を求めることができる。

［ステップＳ１３５］流脈線計算部１３０は、ベクトルｒ_k+1の心臓内部での位置判定を行う。この処理を行うのは、求めたベクトルｒ_k+1は有限の時間幅で誤差を含んでおり、心筋内部に入り込んでしまう場合があるためである。サブルーチンとして位置判定を行った結果、位置判定の結果を示すステータスが返される。例えば心筋を跨いだり、心筋内部に入り込んだりしていなければ正常終了ステータスである「０」または「２」が返される。ベクトルｒ_k+1が解析対象の流体の要素内の位置（例えば心房・心室中）を示している場合、ステータス「０」が返される。ベクトルｒ_k+1が解析対象の流体の要素外の位置（例えば動脈中）を示している場合、ステータス「２」が返される。この処理の詳細は後述する（図１７参照）。

［ステップＳ１３６］流脈線計算部１３０は、判定結果が、正常終了を示すステータス「０」か否かを判断する。ステータスが「０」であれば、処理がステップＳ１４０に進められる。ステータスが「０」でなければ、処理がステップＳ１３７に進められる。

［ステップＳ１３７］流脈線計算部１３０は、判定結果が、終了を示すステータス「２」か否かを判断する。ステータス「２」が返される場合とは、計算の途中で点Ｐ_ijが流体境界を超えて、大動脈のような系外に出た場合である。ステータスが「２」であれば、処理がステップＳ１３８に進められる。ステータスが「２」でなければ、処理がステップＳ１３９に進められる。

［ステップＳ１３８］流脈線計算部１３０は、点Ｐ_ijが系外に出た場合、流脈線上の点Ｐ_ijのｊ＝１からｊまでの領域内判定フラグを「Ｆ」に設定する。領域内判定フラグ「Ｆ」は、点Ｐ_ijが解析領域外に出たことを示している。その後、流脈線計算部１３０は、時間発展計算処理を終了する。このように終了ステータスが「２」の場合、系外に出たとして、領域内判定フラグを「Ｆ」に変えられる。線分を描画するとき、線分は点Ｐ_i1，Ｐ_i2，Ｐ_i3，・・・，，Ｐ_iNの順番に点と点を曲線で結んで描画される。そのため、Ｐ_ijが領域外に出た場合は、粒子発生源から射出したそれ以前の点も描画をする根拠が失われる。そこで流脈線計算部１３０は、ｊ＝１，２，・・・，ｊまでの領域内判定フラグをすべて「Ｆ」に変えて処理を終了している。

［ステップＳ１３９］流脈線計算部１３０は、正常終了しなかった場合（ステータスが「０」または「２」で無い場合）、コントロールパラメータである時間刻み幅を小さくし、さらに計算を複数回に分けて、ステップＳ１３２〜Ｓ１３７の処理を繰り返す。流脈線計算部１３０は、時間刻み幅を小さくしていき、正常に終了したら、処理をステップＳ１４０に進める。

なお、ベクトルｒ_k+1が予測球外に出てしまっているときも、正常終了しなかった場合に相当する。この場合、流脈線計算部１３０は、例えば予測球の半径を大きくして再計算することで、ベクトルｒ_k+1が予測球外に出てしまうことを抑止できる。

［ステップＳ１４０］流脈線計算部１３０は、ベクトルｒ_k+1を保存し、保存した値をｋ＋１番目の繰り返し計算の初期値とする。
［ステップＳ１４１］流脈線計算部１３０は、ステップＳ１３２〜Ｓ１４０の処理が１回終わるごとに、ｋに１を加算し、処理を繰り返す。流脈線計算部１３０は、ｋ＝Ｎ_divの処理が終了すると、処理をステップＳ１４２に進める。

［ステップＳ１４２］流脈線計算部１３０は、最終的に求まったベクトルｒ_Ndivを、ベクトルｒ_k+1として保存する。
図９、図１０に示した処理により、流脈線上の点の座標が更新されていく。

図１１は、流脈線のデータ例を示す図である。図１１に示すように、解析対象の時刻ごとに、各流脈線の点の座標値が設定されている。時刻ｔ＝ｔ₀では、各流脈線の上の点には、初期値としてすべて粒子発生源の座標値が設定されている。時刻ｔ＝ｔ₁において最初に放出された粒子の位置を示す点の座標値が更新される。その後、時刻が進むごとに、新たな粒子が放出され、新たに放出された粒子と既に放出されている粒子それぞれの位置を示す点の座標値が更新される。図１１では、前回の時刻から位置が変更された点の座標値にアンダーラインを示している。

次に、時刻ｔ＝ｔ_kにおける流体部分の構造情報（グリッド座標）、および速度場の情報、弾性体（心筋）部分の構造情報（グリッド座標）とを求める処理（ステップＳ１３２，Ｓ１３３）について詳細に説明する。流体部分の構造情報（グリッド座標）および速度場の情報は、式（３）、式（４）、式（５）、式（６）での計算に使用される。流体の構造情報であるベクトルＢ（ベクトルｒ，ｔ_k）と速度場ベクトルｖ（ベクトルｒ，ｔ_k）である。弾性体の構造情報は、流脈線が心筋と交わっているかの判定に使われ、ベクトルＭ（ベクトルｒ，ｔ_k）である。

データが出力されているのは時刻ｔ＝ｔ_i，ｔ_i+1の場合だけである。そのため、ルンゲ・クッタ法で計算する中間時刻ｔ_kにおいて、グリッド座標、および速度場は定義されていないので、流脈線計算部１３０は、出力ファイルの速度場から場を近似的に求める。このとき留意すべきは、シミュレーション実行時に、時々刻々とグリッドの位置を動かしていくことである。グリッドの位置の決め方には任意性があるが、心臓など境界自体が移動する問題ではＡＬＥ（Arbitrary Lagrangian-Eulerian）法がよく利用される。ＡＬＥ法では、シミュレーションに使用する座標自体が、記述する偏微分方程式の解の精度を落とさないように独立に決定される。決定には一般には偏微分方程式が使用される。しかし、データ解析の立場ではグリッドの位置を決定する支配方程式を知ることはできず、与えられたグリッド点の出力値だけが分かっている。このとき、グリッド点の位置が時間に対して連続的に変化する。

図１２は、グリッド点の位置が時間に対して連続的に変化する様子を示す図である。図１２では、横軸に時刻、縦軸にグリッド点（頂点）のＸ座標値を示している。図１２に示すように、グリッド点は、時間進行に伴い位置が連続的に変化する。そこで流脈線計算部１３０は、グリッド点の位置が時間に対して連続的に変化することを利用して、補間法を用いて任意時刻におけるグリッド位置を推定する。

図１３は、時刻ｔ＝ｔ_kにおける場の情報を求める処理の手順の一例を示すフローチャートである。以下、図１３に示す処理をステップ番号に沿って説明する。
［ステップＳ１５１］流脈線計算部１３０は、時刻ｔ_kにおける場の情報をメモリに設定する。設定される情報は、ベクトルＢ（ベクトルｒ，ｔ_i）、ベクトルｖ（ベクトルｒ，ｔ_i）、ベクトルＢ（ベクトルｒ，ｔ_i+1）、およびベクトルｖ（ベクトルｒ，ｔ_i+1）である。

なお、時刻ｔ_kは、ｔ_i≦ｔ_k≦ｔ_i+1を満たす。そしてベクトルＢ（ベクトルｒ，ｔ_i）、ベクトルｖ（ベクトルｒ，ｔ_i）、ベクトルＭ（ベクトルｒ，ｔ_i）、ベクトルＢ（ベクトルｒ，ｔ_i+1）、ベクトルｖ（ベクトルｒ，ｔ_i+1）、およびベクトルＭ（ベクトルｒ，ｔ_i+1）は既知である。これらのベクトルの値は、ステップＳ１２４とステップＳ１２５において既にファイルから読み込まれている。

流脈線計算部１３０は、はじめに、流体部分の構造を規定するグリッド点について以下の処理を行う。ただし、計算量の削減のため、点Ｐ_i+1jの最大移動距離を理論的に求め、その最大移動距離を半径Ｒとする球内部のグリッド点のみを対象とする。この球を以後予測球と呼ぶ。

［ステップＳ１５２］流脈線計算部１３０は、予測球半径Ｒを決定する。この処理の詳細は後述する（図２０参照）。
［ステップＳ１５３］流脈線計算部１３０は、予測球半径Ｒ内部にある流体のグリッド点を検索し、予測球内部に存在するグリッド点の数をＮ_B,elemに設定する。

［ステップＳ１５４］流脈線計算部１３０は、予測球半径Ｒ内部にあるＮ_B,elem個のグリッド点すべてに対して、ステップＳ１５５〜Ｓ１６０の処理を行う。
［ステップＳ１５５］流脈線計算部１３０は、グリッド点ｉが境界面上か否かを判断する。境界面上であれば、流脈線計算部１３０は、処理をステップＳ１５９に進める。境界面上でなければ、流脈線計算部１３０は、処理をステップＳ１５６に進める。

［ステップＳ１５６］流脈線計算部１３０は、グリッド点が境界面上の点ではないとき、グリッド点ｉの平均移動速度ベクトルｖ_ave,iを計算する。この計算は、時刻ｔ_kにおけるグリッド点の座標を補間して求めることとなる。このとき、各々のグリッド点が出力される時間は短いので、１次式で近似することができる。従って、あるＩＤが与えられたグリッド点の時刻ｔでの座標ベクトルｒ_ID（ｔ）についての平均移動速度ベクトルｖ_aveを以下の式で計算できる。

［ステップＳ１５７］流脈線計算部１３０は、平均移動速度ベクトルｖ_aveを用いて、グリッド点ｉの位置ベクトルｒ_ID（ｔ）を計算する。位置ベクトルｒ_ID（ｔ）は、以下の式で計算できる。

［ステップＳ１５８］流脈線計算部１３０は、非境界面上のグリッド点の速度場ベクトルＶ（ベクトルｒ_ID（ｔ），ｔ）を計算する。この計算は、以下の式（９）に基づいて行われる。

［ステップＳ１５９］流脈線計算部１３０は、境界面上のグリッド位置を示すベクトルｒ_ID（ｔ）を以下の式（１０）に基づいて計算する。

［ステップＳ１６０］流脈線計算部１３０は、境界面上のグリッド位置を示すベクトルｒ_ID（ｔ）の速度場ベクトルＶ（ベクトルｒ_ID（ｔ），ｔ）を以下の式（１１）に基づいて計算する。

［ステップＳ１６１］流脈線計算部１３０は、ステップＳ１５５〜Ｓ１６０の処理を１回行うごとにｉに１を加算し、ステップＳ１５５〜Ｓ１６０の処理を繰り返す。流脈線計算部１３０は、ｉ＝Ｎ_B,elemに対する計算が終了すると、予測球内部の全てのグリッド点における速度場が計算されることになる。

［ステップＳ１６２］流脈線計算部１３０は、計算すべき点ベクトルｒ（ｔ）が含まれる要素番号ＩＤ₀を求める。
［ステップＳ１６３］流脈線計算部１３０は、計算すべき点ベクトルｒ（ｔ）での速度場ベクトルＶ（ベクトルｒ（ｔ），ｔ）を要素番号ＩＤ₀に含まれるグリッド点の情報ｒ_ID0（ｔ）の情報から以下の式（１２）を使って計算する。

その後、時刻ｔ＝ｔ_kにおける場の情報を求める処理を終了する。
図１３に示した処理では、ステップＳ１５７において、非境界面上の位置ベクトルｒ_ID（ｔ）の計算とは別に、ステップＳ１５９において境界面上のグリッドの位置ベクトルｒ_ID（ｔ）を計算している。これは、境界面Ｓ（ｔ）にはナビエ・ストークス方程式に滑りなし境界条件が定められており、拘束条件となるためである。そのため流脈線計算部１３０は、計算対象のグリッド点が境界面Ｓ（ｔ）上の点であるかどうか判定し、境界面Ｓ（ｔ）上の点であれば、境界条件を満たすように決めることとなる。そのために流脈線計算部１３０は、以下の計算をする。

境界面Ｓ（ｔ）上のグリッド点の位置と、その速度場から構成される位相空間上の点を（ベクトルｒ_ID（ｔ），ベクトルｖ_ID（ｔ））とすると、時刻ｔ_iおよびｔ_i+1では値が出力されている。そのため、（ベクトルｒ_ID（ｔ_i），ベクトルｖ_ID（ｔ_i））、（ベクトルｒ_ID（ｔ_i+1），ベクトルｖ_ID（ｔ_i+1））と判明する。いま、粘性流体では、ナビエ・ストークス方程式において滑りなし境界条件を設定するので、グリッド位置ベクトルｒ_ID（ｔ）と速度場ベクトルｖ_ID（ｔ）について以下の関係が成立する。

従って、データは合計４つあるため、境界面Ｓ（ｔ）上のグリッド点については、ベクトルｒ_ID（ｔ）を時間について３次式で定めることにより、これらの条件を満たすようにすることができる。従って式（１０）が得られる。

係数ベクトルａ_iは、ベクトルｒ_ID（ｔ）の成分をξ_ID（ξ＝ｘ，ｙ，ｚ）、その微分の成分をｖξ（ξ＝ｘ，ｙ，ｚ）とすると、

となる。ここで、

とした。
このようにして、任意の時刻における境界面Ｓ（ｔ）上の位置ベクトルｒ_ID（ｔ）の集合を評価できる（図１３のステップＳ１５９）。従って、流体におけるすべてのグリッド点の位置を計算することで流体の構造情報であるベクトルＢ（ベクトルｒ，ｔ_k）を決定することができる。

次に流脈線計算部１３０による、任意の時刻における速度場の計算方法について説明する。速度場はそれぞれ、Ｘ成分、Ｙ成分、Ｚ成分に分けることができる。それぞれの成分についてはスカラー場である。従って、速度場のｊ＝Ｘ，Ｙ，Ｚ成分をＶ_j（ベクトルｒ，ｔ）とし、位置ベクトルｒに十分近く、出力ファイルとしてデータが与えられている時空上の点（ベクトルｒ_{k i}，ｔ_i）の周りで１次の項までテイラー展開し、２次以降の項を無視すると、

となる。ここで、インデックスｋはグリッドＩＤを表し、インデックスｉは出力ファイルが存在するｉ番目の時刻を意味する。∇Ｖ_jは有限要素法における各要素の頂点での速度勾配であるため、

として求めることができる。ここで、Ｎ_kは構造関数と呼ばれ、多項式で与えられる。Ｖ_jkは各頂点における速度場の値であり、ファイルに出力されている。従って、式（２１）における速度場の時間微分項∂Ｖ_j／∂ｔを差分で評価すれば、非境界面における速度場を評価することができる。すなわち、式（９）で速度場を評価することができる（図１３のステップＳ１６０）。場の情報と同様に、心筋の構造情報であるベクトルＭ（ベクトルｒ，ｔ）についても、任意の時刻でのグリッド点情報を計算することができる。

図１４は、時刻ｔ＝ｔ_kにおける弾性体場の情報を求める処理の手順の一例を示すフローチャートである。以下、図１４に示す処理をステップ番号に沿って説明する。
［ステップＳ１７１］流脈線計算部１３０は、読み込まれた心筋の構造情報をメモリに設定する。設定される情報は、ベクトルＢ（ベクトルｒ，ｔ_i）、ベクトルｖ（ベクトルｒ，ｔ_i）、ベクトルＭ（ベクトルｒ，ｔ_i）、ベクトルＢ（ベクトルｒ，ｔ_i+1）、ベクトルｖ（ベクトルｒ，ｔ_i+1）、ベクトルＭ（ベクトルｒ，ｔ_i+1）である。

［ステップＳ１７２］流脈線計算部１３０は、予測球半径Ｒに基づいて、予測球内部の心筋の構造を検索し、心筋のグリッド点を検索する。予測球半径Ｒは、流体部分の処理の時と同じ値が用いられる。流脈線計算部１３０は、検出したグリッド点の数を、グリッド点Ｎ_M,elemに設定する。

［ステップＳ１７３］流脈線計算部１３０は、予測球半径Ｒ内部にあるＮ_M,elem個のグリッド点すべてに対して、ステップＳ１７４〜Ｓ１７８の処理を行う。
［ステップＳ１７４］流脈線計算部１３０は、グリッド点ｉが境界面Ｓ（ｔ）の点であるか否かを判断する。境界面Ｓ（ｔ）上であれば、流脈線計算部１３０は、処理をステップＳ１７７に進める。境界面Ｓ（ｔ）上でなければ、流脈線計算部１３０は、処理をステップＳ１７５に進める。

［ステップＳ１７５］流脈線計算部１３０は、グリッド点が非境界面Ｓ（ｔ）上のとき、個々のグリッド点について流体部分と同様にしてグリッド点ｉの平均移動速度ベクトルｖ_ave,iを式（７）で計算する。

［ステップＳ１７６］流脈線計算部１３０は、グリッド点ｉの位置ベクトルｒ_ID（ｔ）を式（８）に従って計算する。
［ステップＳ１７７］流脈線計算部１３０は、境界面Ｓ（ｔ）上の点であるなら、流体部分とグリッド位置が同じになるので、対応するグリッドを流体部分から探す。

［ステップＳ１７８］流脈線計算部１３０は、見つかったグリッドの座標を境界面Ｓ（ｔ）上の弾性体（心筋）の座標とする。ただし、この処理は流体部分の境界面Ｓ（ｔ）の点が先に計算されることを前提としている。

［ステップＳ１７９］流脈線計算部１３０は、ステップＳ１７４〜Ｓ１７８の処理を１回行うごとにｔに１を加算し、ステップＳ１７４〜Ｓ１７８の処理を繰り返す。流脈線計算部１３０は、ｉ＝Ｎ_M,elemに対する計算が終了すると、時刻ｔ＝ｔ_kにおける弾性体の構造情報を求める処理を終了する。このようにして、心筋の構造情報ベクトル（ベクトルｒ，ｔ）を求めることができる。

次に、心臓内部での位置判定について詳細に説明する。
ここでは、時間発展の結果得られたベクトルｒ_k+1の心臓中の位置を判定する手続きについて説明する。ルンゲ・クッタ法などの方法は有限の誤差を含むため、計算の結果得られた位置ベクトルｒ_k+1が本来であればありえない位置になることがある。

図１５は、計算の結果得られる可能性のある位置の例を示す図である。心臓４０の右心系４１と左心系４２の内部が、解析対象の流体が存在する領域である。右心系４１には、右心房や右心室が含まれる。第２の実施の形態では、左心系４２には、左心房や左心室が含まれる。左心系４２に繋がっている動脈４３内の血流は、解析対象外である。ただし、大動脈の一部をシミュレーションに含めてもよい。

時刻ｔ＝ｔ_kで点Ｐ_kjが流体中に存在するとき、時間発展により移動したときの移動先は、図１５に示すように５通りの分類が可能である。移動した点には、移動先に応じて状態変数（ステータス）が設定される。心筋壁を超えず、解析対象の流体の要素内にある場合の状態変数は「０」である。心筋壁を超え、解析対象の流体の要素外にある場合の状態変数は「１」である。心筋壁を超えず、解析対象の流体の要素外にある場合の状態変数は「２」である。心筋壁を超え、解析対象の流体の要素内にある場合の状態変数は「３」である。心筋壁を超えてはいないが、心筋と解析対象の流体の要素との境界上にある場合の状態変数は「４」である。

図１６は、状態変数の真理値表を示す図である。なお、どちらの場合も移動前の点をＰとし、点Ｐは必ず流体中に存在していると仮定し、移動後の点Ｑの判定のみ行うものとする。どのような状態変数になるのかは、２つの判定を行うことで分類可能である。

１つ目は対象となる点Ｑが流体中に存在するかどうかを判定する流体判定である。点Ｑが流体中に存在すればＴを返し、存在しなければＦを返す。
２つ目は移動元の点Ｐと移動後の点Ｑを結んだ直線ＰＱが心筋または心筋表面と交わるかどうかを判定する直線判定である。直線ＰＱが心筋（表面）と交わればＴを返し、交わらない場合はＦを返す。また直線判定では、同時に交点の個数も返す。

流脈線計算部１３０は、点Ｑが流体中に存在し、直線ＰＱが心筋と交わらない場合は、正常に移動したものとして状態変数「０」を設定する。
点Ｑが流体中に存在せず、直線ＰＱが心筋（表面）と交わる場合、点Ｑは（１）心筋を突き破って系外に出た場合と（２）心筋内部にめり込んだ場合との２通りが存在する。どちらも再計算をすることになるため、流脈線計算部１３０は、状態変数として「１」を設定する。

点Ｑは流体中に存在せず、線分ＰＱが心筋（表面）と交わらない場合は、大動脈などを通してシミュレーション系外に出たものとして、流脈線計算部１３０は状態変数「２」を設定する。

点Ｑが流体中に存在しても、左心房から右心房へ移動するなどありえない移動をする場合では、流体判定はＴが返るが、直線判定もＴであり、交点は必ず複数になる。そのため、交点の数が２以上の場合、流脈線計算部１３０は状態変数「３」を設定する。

流体判定がＴであり、直線判定がＴの場合でも、交点数が１の場合は流体と心筋の境界面に達したとして、流脈線計算部１３０は状態変数「４」を設定する。
このような状態判定の処理手順を以下に示す。

図１７は、状態判定の処理手順の一例を示すフローチャートである。以下、図１７に示す処理をステップ番号に沿って説明する。
［ステップＳ１８１］流脈線計算部１３０は、時間発展後の座標ベクトルｒ_k+1が流体中に存在するかどうか判定する。この処理の詳細は後述する（図１８参照）。

［ステップＳ１８２］流脈線計算部１３０は、移動ベクトルｄｒ＝ベクトルｒ_k+1−ベクトルｒ_kが作る有限の長さの線分が心筋と交点を持つか判定すると共に、この線分が心筋表面と交わる交点の数を計数する。この処理の詳細は後述する（図１９参照）。

［ステップＳ１８３］流脈線計算部１３０は、図１６に示した真理値表に基づいて、状態を判定する。すなわち、ステップＳ１８１，Ｓ１８２の結果は、それぞれ真偽値Ｔ（真、Ｔｒｕｅ）、Ｆ（偽、Ｆａｌｓｅ）の２値で返される。ステップＳ１８３の結果は０以上の整数値となる。流脈線計算部１３０は、真理値表を参照し、ステップＳ１８１〜Ｓ１８３の戻り値に対応する心臓の状態変数として、「０」〜「４」のいずれかの値を位置判定結果とする。

このようにして、位置判定が行われ、状態変数が決定される。
次に、時間発展後の位置が流体か否かの判定処理（ステップＳ１８１）について詳細に説明する。

図１８は、時間発展後の位置が流体内か否かの判定処理の手順の一例を示すフローチャートである。図１８の例では、座標ベクトルｒ_k+1が流体中に存在するか判定する。以下、図１８に示す処理をステップ番号に沿って説明する。

［ステップＳ１９１］流脈線計算部１３０は、時刻ｔ_kでの座標ベクトルｒ_kを取得する。
［ステップＳ１９２］流脈線計算部１３０は、時刻ｔ_k+1での座標ベクトルｒ_k+1を取得する。

［ステップＳ１９３］流脈線計算部１３０は、予測球半径の値Ｒを取得する。
［ステップＳ１９４］流脈線計算部１３０は、座標ベクトルｒ_kを中心とした半径Ｒの予測球内部の要素のリスト（流体要素リストＬ_f）を取得する。

［ステップＳ１９５］流脈線計算部１３０は、流体要素リストＬ_fの要素数をｓｉｚｅ_fとする。
［ステップＳ１９６］流脈線計算部１３０は、結果（ｒｅｓｕｌｔ）の初期値として「Ｆ（Ｆａｌｓｅ）」に設定する。

［ステップＳ１９７］流脈線計算部１３０は、ｉ＝１，２，・・・，ｓｉｚｅ_fまで、要素リスト内の要素Ｌ_iに対してステップＳ１９８の処理を繰り返す。
［ステップＳ１９８］流脈線計算部１３０は、要素リスト内の要素Ｌ_iに対して、要素Ｌ_iが座標ベクトルｒ_k+1を含むかどうかを判定する。座標ベクトルｒ_k+1を含む場合、流脈線計算部１３０は、処理をステップＳ２００に進める。座標ベクトルｒ_k+1を含まない場合、流脈線計算部１３０は、処理をステップＳ１９９に進める。

［ステップＳ１９９］流脈線計算部１３０は、ステップＳ１９８の処理を１回行うごとにｉに１を加算し、ステップＳ１９８の処理を繰り返す。そして流脈線計算部１３０は、ｉ＝ｓｉｚｅ_fの処理が終了したら、時間発展後の位置が流体内か否かの判定処理を終了する。

［ステップＳ２００］流脈線計算部１３０は、要素Ｌ_iの要素番号ＩＤをメモリに保存する。
［ステップＳ２０１］流脈線計算部１３０は、ｒｅｓｕｌｔをＴｒｕｅに変更する。

このようにして、移動先の点の座標がいずれかの要素に含まれていた場合、座標ベクトルｒ_k+1は流体中に存在するとして戻り値がＴ（Ｔｒｕｅ）に設定される。しかし、含まれていない場合は戻り値がＦ（Ｆａｌｓｅ）に設定される。

次に、予測球内部にある移動ベクトルｄｒが通過する弾性体要素の検索処理（ステップＳ１８２）について詳細に説明する。
図１９は、移動ベクトルｄｒが通過する弾性体要素の検索処理の手順の一例を示すフローチャートである。以下、図１９に示す処理をステップ番号に沿って説明する。

［ステップＳ２１１］流脈線計算部１３０は、時刻ｔ_kでの座標ベクトルｒ_kを取得する。
［ステップＳ２１２］流脈線計算部１３０は、時刻ｔ_k+1での座標ベクトルｒ_k+1を取得する。

［ステップＳ２１３］流脈線計算部１３０は、取得した情報から移動ベクトルｄｒ＝ベクトルｒ_k+1−ベクトルｒ_kを計算する。これにより、点の移動経路を示す線分が定義される。

［ステップＳ２１４］流脈線計算部１３０は、予測球半径Ｒの値を取得する。
［ステップＳ２１５］流脈線計算部１３０は、予測球半径内部にある弾性体の要素リストＬ_eを作成する。

［ステップＳ２１６］流脈線計算部１３０は、流体要素Ｌ_eの要素数をｓｉｚｅ_eとする。
［ステップＳ２１７］流脈線計算部１３０は、結果（ｒｅｓｕｌｔ）の初期値として「Ｆ（Ｆａｌｓｅ）」に設定する。

［ステップＳ２１８］流脈線計算部１３０は、ｉ＝１，２，・・・，ｓｉｚｅ_eまで、要素リスト内の要素Ｌ_eに対してステップＳ２１９の処理を繰り返す。
［ステップＳ２１９］流脈線計算部１３０は、ｉ番目の要素Ｌ_iが線分ベクトルｄｒと交わるか否かを判断する。要素Ｌ_iは多面体であるため、流脈線計算部１３０は、要素Ｌ_iに含まれるすべての面に対して線分ベクトルｄｒと交点を求める。要素内のすべての面に対して、交点が１点も存在しなければ、流脈線計算部１３０は、その要素とは交わらないと判定する。交わらない場合、流脈線計算部１３０は、処理をステップＳ２２０に進める。交わる場合、流脈線計算部１３０は、処理をステップＳ２２１に進める。

［ステップＳ２２０］流脈線計算部１３０は、ステップＳ２１９の処理を１回行うごとにｉに１を加算し、ステップＳ２１９の処理を繰り返す。そして流脈線計算部１３０は、ｉ＝ｓｉｚｅ_eの処理が終了したら、移動ベクトルｄｒが通過する弾性体要素の検索処理を終了する。

［ステップＳ２２１］流脈線計算部１３０は、要素Ｌ_iの要素番号ＩＤをメモリに保存する。
［ステップＳ２２２］流脈線計算部１３０は、線分を示すベクトルｄｒと心筋表面との交点数Ｚ（Ｚは１以上の整数）を求め、メモリ１０２に保存する。

［ステップＳ２２３］流脈線計算部１３０は、ｒｅｓｕｌｔをＴｒｕｅに変更する。
このようにして、交点が１点も存在しない場合、交わらないとしてＦが返され、交点が１点以上存在する場合、交わるとしてＴが返される。交わる場合は交点を持つ要素の番号と、心筋表面との交点数Ｚが保存される。

以上の処理により、点の移動先の状態を適切に判定することができる。この際、移動先の判定対象として予測球内の要素を対象とすることで、処理を効率化することができる。以下、予測球半径Ｒの決定方法について詳細に説明する。

流脈線計算部１３０は、予測球は時間発展前の座標ベクトルｒ_kが時間刻み幅Δｔの間にどれだけの距離を動きえるかに基づいて決定する。４次のルンゲ・クッタ法の場合、式（２）式から次の不等式が成立する。

従って、半径Ｒを

と定義すれば、半径Ｒは流脈線上の点Ｐが時間Δｔに進む最大距離を表す。そして、時間Δｔ後にも点Ｐは必ずこの球内部に存在する。また、中間ベクトルｖ_iも同様にして、半径Ｒの球内部の点であることを次のようにして示すことができる。すなわち式（４）と式（５）において、

が成立する。そのため、中間ベクトルが指し示す座標も半径Ｒの球内部の点になるのである。同様のことを式（３）、式（６）について行えば、すべての中間ベクトルｖ_iも同様にして、半径Ｒの球内部の点であることが分かる。このような半径Ｒを予測球半径とする。流脈線計算部１３０は、予測球半径を設定するため流体の速度場から最大値を求める。単純には、シミュレーションを行った領域全部のグリッド点の中から最大値を求めればよいが、シミュレーション対象の領域を幾つかの適切な領域に分割してから、局地的な情報から最大値を求めることもできる。

ナビエ・ストークス方程式の性質から、速度の最大値は境界面Ｓ（ｔ）上にはありえないので、速度の最大値は中間時刻でも式（９）で補間される。式（９）は時間に関して２次式であるため、データが出力されている二つの時刻ｔ＝ｔ_i、ｔ＝ｔ_i+1のデータから、その２点を通る放物線の最大値から速度の最大値を求めることが可能である。近似的には、式（９）において２次の項はテイラー展開の性質上、２次の微小量になる。そのため、１次の項のみを考え、ｔ＝ｔ_i、ｔ＝ｔ_i+1におけるそれぞれの時刻における速度場の最大値を比較し、大きい方を速度場の最大値としてもよい。従って流脈線計算部１３０は、求めた速度場のノルム最大の点を最大速度｜ベクトルｖ_max｜として設定する。その後、ルンゲ・クッタ法の時間刻み幅Δｔを設定すると、最大移動距離はＲ＝｜Δｔｖ_max｜として定まる。

しかし、このように定めた予測球半径Ｒは時間刻み幅Δｔが小さくなりすぎると、グリッドの幅（隣接するグリッド点間の距離）より小さくなってしまう。そのため、グリッド点の離散化に伴う最小値Ｒ_minが存在する。例えば、このＲ_minは経験的な最小値として初期値をＲ_min＝０．００１［ｍ］に設定する。グリッドの統計分析から最小値を求めることも可能である。

予測球の半径Ｒが、このようにして求めたＲ_minより小さくなった場合、予測球の半径はＲ＝Ｒ_minとして、予測球内部に要素が一つも存在しない状況を回避することができる。以上のような予測球半径Ｒの決定処理を整理すると、図２０のような手順となる。

図２０は、予測球半径の決定処理の手順の一例を示すフローチャートである。以下、図２０に示す処理をステップ番号に沿って説明する。
［ステップＳ２３１］流脈線計算部１３０は、要素の辺の長さと速度場の統計解析から予測球半径の最小値Ｒ_minを決定する。

［ステップＳ２３２］流脈線計算部１３０は、時刻ｔ_iでの速度場ベクトルｖ_iの値を取得する。
［ステップＳ２３３］流脈線計算部１３０は、時刻ｔ_iでの速度の最大値｜ベクトルｖ_i,max｜を取得する。

［ステップＳ２３４］流脈線計算部１３０は、時刻ｔ_i+1での速度の最大値｜ベクトルｖ_i+1,max｜を取得する。
［ステップＳ２３５］流脈線計算部１３０は、｜ベクトルｖ_i,max｜が｜ベクトルｖ_i+1,max｜以上かどうかを判断する。流脈線計算部１３０は、｜ベクトルｖ_i,max｜が｜ベクトルｖ_i+1,max｜以上であれば、処理をステップＳ２３６に進める。流脈線計算部１３０は、｜ベクトルｖ_i,max｜が｜ベクトルｖ_i+1,max｜未満であれば、処理をステップＳ２３７に進める。

［ステップＳ２３６］流脈線計算部１３０は、｜ベクトルｖ_i,max｜を｜ベクトルｖ_max｜に設定する。流脈線計算部１３０は、その後、処理をステップＳ２３８に進める。
［ステップＳ２３７］流脈線計算部１３０は、｜ベクトルｖ_i+1,max｜を｜ベクトルｖ_max｜に設定する。流脈線計算部１３０は、その後、処理をステップＳ２３８に進める。

［ステップＳ２３８］流脈線計算部１３０は、ルンゲ・クッタ法の時間刻み幅ｄｔを取得する。
［ステップＳ２３９］流脈線計算部１３０は、｜ベクトルｖ_max｜ｄｔを計算し、計算結果を予測球半径Ｒに設定する。

［ステップＳ２４０］流脈線計算部１３０は、予測球半径Ｒが、予測球半径の最小値Ｒ_minより小さいか否かを判断する。流脈線計算部１３０は、予測球半径Ｒが最小値Ｒ_minより小さい場合、処理をステップＳ２４１に進める。流脈線計算部１３０は、予測球半径Ｒが最小値Ｒ_min以上の場合、予測球半径決定処理を終了する。

［ステップＳ２４１］流脈線計算部１３０は、最小値Ｒ_minを予測球半径Ｒに設定する。その後、流脈線計算部１３０は、予測球半径決定処理を終了する。
以上のような処理により、適切な予測球半径Ｒを決定することができる。

なお、第２の実施の形態では、以下のような処理によって、処理の高速化を図ることができる。
［分割法による計算の高精度化と高速化］
図９のステップＳ１２７では、出力ファイルが与えられた時間ｔ_i≦ｔ≦ｔ_i+1の時間分割をさらにＮ_div分割している。この分割の仕方について詳細に説明する。

流脈線の計算では一回の流脈線上の点が関与する要素は全体のごく一部である。そのため、記憶容量削減と計算時間の短縮のため予測球を用いる。従って、流脈線の計算コストは予測球半径Ｒ³に比例して増大する。このことは、以下のように説明できる。

計算対象となる要素の数Ｎ_elemは、空間的な要素数の密度を座標ベクトルｒの関数としてρ（ベクトルｒ）とすると、以下の式で与えられる。

密度ρ（ベクトルｒ）がおおよそ一様であると仮定すれば、ρ（ベクトルｒ）＝ρ₀となる。そのため、Ｎ_elem∝Ｒ³となる。
ここで時間ｔ_i≦ｔ≦ｔ_i+1をＮ_div分割することを考える。Ｎ_div分割すると、ルンゲ・クッタ法による時間発展回数はＮ_div回となる。一方で、予測球半径は式（２４）から決定される。そのため、時間刻み幅がＮ_div分の１になると、予測球半径もＮ_div分の１になる。１回の計算量は予測球半径Ｒ³に比例するため、Ｎ_div ^-3になる。この計算をＮ_div回繰り返すため、合計の計算量はＮ_div ^-2になる。この計算の概要を模式化したのが図２１である。

図２１は、時間分割による計算量削減の概念を示す図である。図２１では、点Ｐの移動先となる点Ｑを計算する場合を想定している。時間分割せずに時間刻み幅を長くすると、点Ｐを中心とする大きな予測球４５内を解析することとなるが、時間分割を行い時間刻み幅を短くすると、点Ｐにある粒子の移動可能範囲が狭まるため、予測球４５よりも小さな予測球４６内を解析すればよい。すなわち、大きな予測球４５に対して１回計算をするより、時間をＮ_div分割して、経路を細分化したほうが予測球４６が小さくなる。その結果、計算対象になる要素数を減らせ、計算量を削減できることが分かる。

図２２は、時間分割数に応じた合計計算時間の変化の一例を示す図である。図２２のグラフは、横軸が分割数Ｎ_divであり、縦軸が計算時間である。理論上の計算時間５１の遷移を破線で示し、実際の計算時間５２の遷移を実線で示している。実際の計算時間５２は、時間ｔ_i≦ｔ≦ｔ_i+1のトータルの計算時間を測定したものである。Ｎ_div＝０の時は、予測球を使用せず、そのまま全部の要素を計算対象とした場合にかかった時間である。

理論上は分割数Ｎ_divが多くなるほど計算時間は短くなる。しかし、実際には、分割数を多くしすぎると、予測球内の要素のリスト作成など、分割後の単位時間当たり１回ずつ行われる処理の回数が増えてしまい、処理時間が増加に転じる。そのため、分割数Ｎ_divには最適値が存在する。分割数Ｎ_divの最適値は、例えば流脈線の計算を始める前にＮ_divを変えながら、対象となる系に対して予め測定しておく。

［統計分析による予測球半径の最小値の決定法］
流体シミュレーションは有限個の離散点情報で行われる。従って、予測球半径を小さくしすぎると予測球内部に離散点情報が含まれない状況が生じ得る。これは予測球半径に下限値が存在することを意味する。同時に、予測球半径は式（２４）から決定されるため、時間刻み幅にも下限値が存在する。予測球を用いた計算では、予測球半径の下限値Ｒ_minを設定し、その予測球半径が下限値以下の場合になるときは、Ｒ_minを予測球半径として設定する。なお、予測球の下限値以下になった時に、それより大きな値であるＲ_minを使用するので、流脈線上の点が予測球外部に出てしまうことはない。このように、計算を安定に進めるためにＲ_minの決定は重要である。また、Ｒ_minの値は時間分割数Ｎ_divの決定にも関与する。Ｒ_minの値が定まれば、Ｒ₀＝｜Δｔベクトルｖ_max｜が最大移動距離であるため、天井関数を使って以下の式で分割数Ｎ_divを決定できる。

分割数Ｎ_divの値は計算速度と計算精度に関わるため、Ｒ_minの決定は計算速度と計算精度の面でも重要である。
しかし、下限値Ｒ_minの決め方には工夫がいる。具体的には、確率モデルを導入し、流脈線計算部１３０は、計算実行時に失敗することを許容した投機的な計算を行う。計算に失敗したときには計算が確実に成功するパラメータで計算を実行し、その計算にかかる時間がペナルティとしてかかると考える。そして、ペナルティまで含めて統計的に計算時間を最小にするように下限値Ｒ_minを決定する。詳細を以下に記す。

流脈線計算部１３０は、初めにペナルティとして用いる最悪値としての半径Ｒ_wを決める。これはシミュレーションに用いられた要素の辺の長さの最長のものを設定する。
図２３は、要素の辺の長さの分布の一例を示す図である。図２３のグラフは、横軸が要素の辺の長さであり、縦軸が、辺の長さごとの該当する辺の存在確率を示す確率密度である。グラフ中には、２種類の心臓シミュレーション＃１，＃２における確率密度関数６１，６２を破線で示している。

図２３の例では、シミュレーションに使用された要素の一部が非常に粗く、最長で０．００８［ｍ］にも達する。そのため、そのまま計算に用いてしまうと分割数Ｎ_divも小さくなり、速度が遅くなる。しかし、殆どの計算は最長の長さでなくても計算できる。

最悪値としての半径Ｒ_wの決定後、流脈線計算部１３０は、計算に使用する下限値Ｒ_minを決定する。この際、流脈線計算部１３０は、統計的な計算コストを、計算の失敗を前提として計算する。そして、流脈線計算部１３０は、投機的な実行を考慮に入れ、失敗した時のペナルティまで含めて計算量が最小になるような半径を、下限値Ｒ_minとして選択する。これは、仮定する確率モデルに極端に依存する話であるが、簡単な例を用いて説明する。まず、流脈線上の点は、すべての要素を等確率で解析対象になると仮定する。すると、図２３の確率密度関数６１，６２に従って、計算が成功する。

今、半径Ｒで計算時間がＴだけかかるとする。この計算量はＲ³に比例する。また半径をβＲ（０＜β≦１）に縮小することで、計算はｐ［％］の確率で成功するものとする。計算が成功する場合とは、点の移動先が半径βＲの予測球内にある場合である。計算に失敗した場合、半径Ｒで再計算する。この場合、ペナルティまで含めた計算時間Ｔ’は、

となる。
図２４は、予測球半径に応じた計算コストの変化を示す図である。図２４の縦軸が計算コスト（Ｔ’／Ｔ）であり、横軸が予測球半径Ｒである。図２４に示すように、Ｒ＝０．００３［ｍ］付近でペナルティまで含めたコストが最悪値半径Ｒ_max＝０．００８［ｍ］と比べて約１０％まで落ちることがわかる。そして、ペナルティまで含めて統計的に最小の計算量で済む半径が存在することがわかる。実際には、有限回の計算であるため、すべての要素を等確率に解析対象とするという状態にはならない。そのため、これよりも半径の値は小さくなる。

この場合、最適な予測球半径Ｒが存在することは示せるが、最適値自体は過大評価される傾向がある。そこで、過大評価を抑えることを目的とし、各離散点における速度場のノルムを求め、これに時間刻み幅をかけて移動距離の分布を出すと図２５のようになる。

図２５は、移動距離の確率分布を示す図である。図２５では、横軸が移動距離、縦軸が発生数を示している。同じ形状、同じ大きさの４面体があったとしても、４面体上での速度場の大きさが２倍違えば、対象となる４面体を通過する時間は半分になる。そこで、移動距離の分布から式（２８）を用いて予測球半径Ｒの関数としてペナルティまで含めた計算効率性を図示したのが図２６である。

図２６は、確率分布を仮定した時の計算効率化曲線を示す図である。図２６では、横軸が予測球半径Ｒ、縦軸が計算効率である。計算効率は、値が小さいほど効率がよいことを示す。図２６から計算効率が最もよい予測球半径Ｒを求めると、０．００１２［ｍ］〜０．００１６［ｍ］程度になる。これは経験的な最良値Ｒ_min＝０．００１５±０．０００５［ｍ］のモデル論的な裏付けである。

以上のようにして、予測球半径Ｒの最良値を求めることができる。すなわち、離散点が予測球の外部に出てしまう場合の再計算のペナルティを考慮した、最も効率的な予測球半径Ｒを設定することができる。その結果、流脈線の計算効率が向上する。

［具体的な計算例］
以下、心臓シミュレータによる心筋と冠循環の解析を用いて、流脈線を実際に計算した場合の計算速度について具体的に説明する。

初めに計算に用いたデータについて説明する。シミュレーション結果は心臓の拍動１周期分が出力されており、０．０１［ｓｅｃ］ごとに１００個の心臓の状態のデータが出力されている。心臓は拍動するため、速度場は時間依存して変化するため非定常流である。また、心筋も拍動のため拡張と収縮を繰り返す。従って、心筋面も動き、移動境界問題になる。この系において心臓中の血流の輸送を記述するため、心臓中の複数の点に粒子発生装置を配置し、流脈線を計算する。

可視化装置１００に外部から心筋と流体の情報を読み込ませる。この際、可視化装置１００に、流脈線の本数Ｍと計算回数Ｎを入力する。また可視化装置１００に対して、心臓中の左心室と右心房を流脈線の発生源とするように、発生源の位置を入力する。発生源の配置指示に従い、例えば可視化装置１００は、半径０．０５［ｍ］の球内部の流体部分からランダムに２０個ずつ、流脈線の発生源を配置したものとする。また、計算回数は０．３拍分の３０回とする。なお、時間発展計算は４次のオーダのルンゲ・クッタ法で行った。

以上のような初期条件を設定し、可視化装置１００が流脈線計算を開始する。開始後、ステップ毎に粒子発生源から粒子が発生する。それぞれの粒子の時間Δｔ＝０．０１［ｓｅｃ］後の位置は図８に示した流脈線計算のフローチャートに従って、計算される。計算された流脈線は、モニタ２１に表示される。

図２７は、流脈線の表示例を示す図である。図２７の例では、心臓７０の断面の画像の上に重ねて、複数の流脈線７１が表示されている。
流脈線の計算では、計算量の削減のため、および計算精度の向上のため、予測球が使用されている。予測球を使った計算では、時間刻み幅Δｔ＝０．０１［ｓｅｃ］をさらにＮ_div分割して細分化した。Ｎ_divの値は、統計的に計算量が最小になるように自動的に最適値が決定される。このシミュレーション結果の場合、Ｎ_divの値はＮ_div＝３〜７程度に設定された。

Ｎ_divの最適化による計算速度を定量的に測定した結果は、図２２に示した通りである。図２２の例では、予測球を使用した場合の効果を示すため、予測球を使用しない場合の流脈線の計算時間も計測している。Ｎ_div＝０のときが、予測球を使用しない場合、すなわちそのまま全部の要素を計算対象とした場合にかかった時間である。

予測球を使った場合（Ｎ_div＝１）の時、つまり、実質上時間分割をしない場合は予測球の構築などに時間がかかる。そのため、予測球を使用しない場合よりも時間がかかる（約１．４７倍）。しかし、Ｎ_div＝２とすると、予測球の構築にかかる固定コストよりも、Ｎ_div ^-2に比例する計算コスト減少のコストが効いてくる。そのため、計算時間は予測球を使用しない場合と比べても約２２．９％まで減少する。Ｎ_divを増やしていくと、計算時間は短くなっていくが、予測球を構築するための固定コストがあるため、性能向上に上限値が存在することが分かる。図２２の例では、Ｎ_div＝４が最良である。そして、分割数Ｎ_divを増やして行くと、固定コストの効果が上回り、計算時間は増大に転じる。Ｎ_div＝１０では最良のケースと比べて性能が約５３％程度まで悪化している。

このように分割数Ｎ_divには最適値が存在するが、実際の計算では流脈線の計算を始める前に最良の分割数が決定されている。また、図２２の計算では、グリッド点数が高々５０，０００点程度であるため、予測球を使わなかった場合の方が、予測球を使いＮ_div＝１の場合よりも計算時間は短い。しかし、シミュレーション系が大規模になれば、予測球を使わない場合は計算量がグリッド点に比例して増えていく。例えば、グリッド数が１０倍になれば、時間も１０倍かかる。しかし、予測球を使用した場合、予測球半径Ｒは速度の最大値と時間刻みだけから決まり、常に一定になるため、計算量も一定になる。このため、シミュレーションを行う系が大きくなるほど、予測球を使用する効果は大きくなる。

さらに、予測球による時間分割のもう一つの効果について述べる。それは精度向上である。４次のルンゲ・クッタ法の場合、式（２）〜式（６）で与えられるが、誤差のオーダはＯ（Δｔ⁴）で与えられることが知られている。そのため、時間をＮ_div分割すると、１回のルンゲ・クッタ演算での誤差はＮ_div ^-4倍になる。時間ｔ_i≦ｔ≦ｔ_i+1まで合計Ｎ_div回のルンゲ・クッタ法の計算をして、誤差がＮ_div倍だけ累積しても、合計の誤差はＮ_div ^-3に留まる。Ｎ_div＝４ならば、１／６４倍になる。

図２８は、時間刻み幅を変えたときの軌跡の精度の変化を示す図である。図２８では、横軸がシミュレーション上の時間、縦軸が離散点のｘ座標である。図２８の例では、時間刻み幅Δｔを「０．０１」、「０．０５」、「０．００１」、「０．０００１」としたときの、特定の離散点のｘ座標の時間変化（軌跡）を表している。

積分後の座標変化量が小さい場合、もしくは積分回数が少ない場合においては、時間刻み幅が大きい場合でも軌跡はほとんど変わらない。しかし、積分回数が大きくなった領域（時刻ｔ＞２．５［ｓｅｃ］）では時間刻み幅Δｔ＝０．０１［ｓｅｃ］の場合の軌跡は、Δｔ＝０．００１［ｓｅｃ］，０．０００１［ｓｅｃ］の軌跡と比べて大幅にずれてしまっている。反対にΔｔ＝０．００１［ｓｅｃ］の軌跡とΔｔ＝０．０００１［ｓｅｃ］の軌跡はよい一致を見せている。このことから、時間刻み幅を小さくすることで、計算誤差の累積を防ぐことができることがわかる。

このように、分割数Ｎ_divの決定は計算速度の面でも、計算精度の面でも重要である。そして、分割数Ｎ_divは予測球半径の最小値Ｒ_minに依存して決まる。実際の計算では図２５のような速度場と要素の辺の長さの両方を考慮した確率モデルを仮定し、確率分布を計算する。そして、その確率分布から式（２８）を用いて計算効率を求め、最良の計算効率を与えるＲ_minを求める。これは図２６に示すように、最小値が存在し、その値はおおよそＲ_min＝０．００１５±０．０００５［ｍ］程度として与えられる。そして、常に最小の計算コストで計算を実行することができる。

〔その他の実施の形態〕
流脈線上の各点はお互いに相互作用せず、独立である。そのため、並列化に適している。そこでＭＰＩ（Message Passing Interface）やＯｐｅｎＭＰ（Multi-Processing）を用いて流脈線の計算を並列化させてもよい。これにより計算速度を向上させることができる。

また第２の実施の形態では、心臓シミュレーションの結果を用いて、血流の流脈線を可視化する例を示したが、他の流体シミュレーションの結果に対しても適用できる。例えば自動車のリアに可変ウィング機構がついている場合に、可変ウィング稼働時の周囲の空気の流れをシミュレーションによって解析することも考えられる。第２の実施の形態に示す可視化装置１００を用いれば、可変ウィング稼働時の周囲の空気の流れを流脈線で可視化することができる。また、可変翼を有する航空機における翼の変形過程での周囲の空気の流れを示す流体シミュレーション結果についても、同様に適用できる。

１ 解析空間
２ 構造体
３ａ，３ｂ 流体領域
４ 解析対象領域
５ 第１の位置
６ 第２の位置
７ 表示画面
８ 粒子発生源
９ａ，９ｂ 流脈線
１０ 流脈線可視化装置
１１ 記憶部
１２ 処理部

Claims (9)

1. 流体シミュレーション上の時間内の複数の解析時刻における複数の粒子の連なりを示す流脈線を計算し、前記流脈線を表示する流脈線可視化装置であって、
   解析空間内の構造体の形状の時間変化を示す構造体情報と、前記解析空間内の流体が存在する領域内の複数の点における、前記流体の速度の空間変化又は／及び時間変化を示す流体情報とを記憶する記憶部と、
   第１解析時刻における第１流脈線に基づいて第２解析時刻における第２流脈線を計算する場合、前記第１流脈線上の第１の位置にある離散点を含む、前記解析空間内の部分領域を、前記離散点についての解析対象領域に決定し、前記流体情報に示される前記解析対象領域内の前記流体の速度に基づいて、前記第２解析時刻における前記離散点上の粒子の移動先を示す第２の位置を計算し、前記構造体情報に示される前記解析対象領域内に存在する前記構造体の情報に基づいて、前記第２解析時刻における前記解析対象領域内の前記構造体の占有領域を特定し、前記第１の位置と前記第２の位置とに基づいて、前記第２流脈線の前記占有領域への侵入の有無を判定し、前記第２流脈線が前記占有領域に侵入しないと判定したとき、前記第２の位置を通過する前記第２流脈線を表示する処理部と、
   を有する流脈線可視化装置。
2. 前記処理部は、さらに、
   前記第２流脈線が前記占有領域に侵入していると判定した場合、前記第１解析時刻と前記第２解析時刻との間の時間帯内に１つ以上の第３解析時刻を設定し、前記第３解析時刻での前記離散点上の粒子の移動先を示す第３の位置を計算し、前記第３の位置に基づいて前記第２の位置を再計算する、
   請求項１記載の流脈線可視化装置。
3. 前記処理部は、前記解析対象領域の決定において、前記第１解析時刻と前記第２解析時刻とを含む時間帯における前記流体の最大速度を求め、前記第１解析時刻と前記第２解析時刻との差、および前記最大速度に基づいて、球状の領域である前記解析対象領域の半径を算出する、
   請求項１または２記載の流脈線可視化装置。
4. 前記処理部は、前記解析対象領域の決定において、前記流体が存在する領域内の前記複数の点の間隔に基づいて、前記解析対象領域の半径の最小値を設定し、算出した半径が前記最小値より小さい場合、前記最小値を前記解析対象領域の半径に決定する、
   請求項３記載の流脈線可視化装置。
5. 前記処理部は、前記解析対象領域の決定において、決定した半径による前記第２の位置の計算が失敗した場合、代わりに半径の最大値を前記解析対象領域の半径として前記第２の位置を再計算すると共に、失敗した場合の前記第２の位置の計算の量も含めて計算量を最小化する半径を算出し、算出した半径が前記最小値より小さい場合、前記最小値を前記解析対象領域の半径に決定する、
   請求項４記載の流脈線可視化装置。
6. 前記処理部は、
   前記解析対象領域の決定において、前記解析対象領域の半径が、前記第１解析時刻と前記第２解析時刻との差を前記最大速度に乗算した値よりも小さいことを許容し、
   前記第２の位置の計算において、前記離散点上の粒子の移動先が前記解析対象領域外となった場合、前記解析対象領域を拡大して前記第２の位置を再計算する、
   請求項３ないし５のいずれかに記載の流脈線可視化装置。
7. 前記処理部は、侵入の有無の判定において、前記第２の位置が前記流体内にあるか否かの第１の判定と、前記第１の位置と前記第２の位置とを結ぶ線分が前記構造体の表面と交わるか否かの第２の判定とを行い、前記第１の判定と前記第２の判定との結果に基づいて、前記第２流脈線の前記占有領域への侵入の有無を判定する、
   請求項１ないし６のいずれかに記載の流脈線可視化装置。
8. 流体シミュレーション上の時間内の複数の解析時刻における複数の粒子の連なりを示す流脈線を計算し、前記流脈線を表示するための流脈線可視化方法であって、
   コンピュータが、
   第１解析時刻における第１流脈線に基づいて第２解析時刻における第２流脈線を計算する場合、前記第１流脈線上の第１の位置にある離散点を含む、解析空間内の部分領域を、前記離散点についての解析対象領域に決定し、
   前記解析空間内の流体が存在する領域内の複数の点における、前記流体の速度の空間変化又は／及び時間変化を示す流体情報を参照し、前記解析対象領域内の前記流体の速度に基づいて、前記第２解析時刻における前記離散点上の粒子の移動先を示す第２の位置を計算し、
   前記解析空間内の構造体の形状の時間変化を示す構造体情報を参照し、前記解析対象領域内に存在する前記構造体の情報に基づいて、前記第２解析時刻における前記解析対象領域内の前記構造体の占有領域を特定し、
   前記第１の位置と前記第２の位置とに基づいて、前記第２流脈線の前記占有領域への侵入の有無を判定し、
   前記第２流脈線が前記占有領域に侵入しないと判定したとき、前記第２の位置を通過する前記第２流脈線を表示する、
   流脈線可視化方法。
9. 流体シミュレーション上の時間内の複数の解析時刻における複数の粒子の連なりを示す流脈線を計算し、前記流脈線を表示する処理をコンピュータに実行させるための流脈線可視化プログラムであって、
   前記コンピュータに、
   第１解析時刻における第１流脈線に基づいて第２解析時刻における第２流脈線を計算する場合、前記第１流脈線上の第１の位置にある離散点を含む、解析空間内の部分領域を、前記離散点についての解析対象領域に決定し、
   前記解析空間内の流体が存在する領域内の複数の点における、前記流体の速度の空間変化又は／及び時間変化を示す流体情報を参照し、前記解析対象領域内の前記流体の速度に基づいて、前記第２解析時刻における前記離散点上の粒子の移動先を示す第２の位置を計算し、
   前記解析空間内の構造体の形状の時間変化を示す構造体情報を参照し、前記解析対象領域内に存在する前記構造体の情報に基づいて、前記第２解析時刻における前記解析対象領域内の前記構造体の占有領域を特定し、
   前記第１の位置と前記第２の位置とに基づいて、前記第２流脈線の前記占有領域への侵入の有無を判定し、
   前記第２流脈線が前記占有領域に侵入しないと判定したとき、前記第２の位置を通過する前記第２流脈線を表示する、
   処理を実行させる流脈線可視化プログラム。

Images