JP2018180707A - 流脈線可視化装置、流脈線可視化方法、および流脈線可視化プログラム - 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 William H. Press, Brian P. Flannery, Saul A. Teukolsky and William T. Vetterling, "Numerical Recipes in C: The Art of Scientific Computing", Cambridge University Press, 1992/10/30, pp.113-117. Joseph E. Flaherty, "Finite Element Analysis Chapter 4 Finite Element Approximation", 2005年4月1日、［２０１７年２月２７日検索］、インターネット<URL:http://www.cs.rpi.edu/~flaherje/pdf/fea4.pdf>

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

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

１つの案では、記憶部と処理部とを有する流脈線可視化装置が提供される。
記憶部は、流体シミュレーションの解析時間の進行に伴い解析空間内を加速度運動する複数のグリッド点の第１時間間隔の複数の第１時刻の各々における位置を示す位置情報と、複数のグリッド点上の流体の複数の第１時刻の各々における速度を示す速度情報とを含む流体情報を記憶する。処理部は、まず流体情報に基づいて、位置情報に表される複数のグリッド点の加速度運動に起因する誤差を補正する補正値を含む計算式により、複数のグリッド点上の流体の複数の第１時刻の各々における速度の時間微分値をそれぞれ算出する。次に処理部は、複数のグリッド点上の流体の複数の第１時刻の各々における速度と速度の時間微分値とに基づいて、解析時間の進行に伴い粒子発生源から順番に出力される複数の粒子の、第１時間間隔より短い第２時間間隔の複数の第２時刻の各々における位置をそれぞれ算出する。そして処理部は、複数の粒子の連なりを示す流脈線の表示情報を生成する。

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

第１の実施の形態に係る流脈線可視化装置の一例を示す図である。 第２の実施の形態のシステム構成例を示す図である。 可視化装置のハードウェアの一構成例を示す図である。 流脈線の計算例を示す図である。 可視化装置の機能を示すブロック図である。 弾性体情報ファイル群の一例を示す図である。 流体情報ファイル群の一例を示す図である。 事前解析ファイルの一例を示す図である。 流脈線可視化処理の手順の一例を示すフローチャートである。 事前解析処理の手順の一例を示すフローチャートである。 ＡＬＥ制御されたグリッド点自身の運動による見かけの力の原理を示す図である。 時間発展計算処理の手順を示すフローチャートの前半部分である。 時間発展計算処理の手順を示すフローチャートの後半部分である。 ルンゲ・クッタ係数表の一例を示す図である。 流脈線のデータ例を示す図である。 グリッド点の位置が時間に対して連続的に変化する様子を示す図である。 時刻ｔ＝ｔ_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」というファイル名が付与される。ファイル名中の「Ｙ」は、シミュレーションの時刻が古い順に付与された昇順の番号である。

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

図８は、事前解析ファイルの一例を示す図である。事前解析ファイル１２１には、すべての時刻における、すべてのグリッド点の補間多項式（後述の式（２））を再現できるように、グリッド番号、時刻、位置座標、速度場（血流のみ）、補間多項式係数が記録される。例えば事前解析ファイル１２１には、出力時刻インデックステーブル１２１ａ、心筋側補間多項式係数テーブル１２１ｂ、および流体側補間多項式係数テーブル１２１ｃが含まれる。

出力時刻インデックステーブル１２１ａには、時刻のインデックス（Time Index）に対応付けて、流体解析のシミュレーション上での時刻（Time）が設定されている。
心筋側補間多項式係数テーブル１２１ｂには、グリッドのＩＤ（GRID ID）と時刻のインデックス（Time Index）と軸のインデックス（Direction Index）の組に対応付けて、軸方向の位置（Position）、心筋側補間多項式に用いる各係数（係数ａ、係数ｂ、係数ｃ、係数ｄ）の値が設定されている。

流体側補間多項式係数テーブル１２１ｃには、グリッドのＩＤ（GRID ID）と時刻のインデックス（Time Index）と軸のインデックス（Direction Index）の組に対応付けて、軸方向の位置（Position）、軸方向の速さ（Velocity）、流体側補間多項式に用いる各係数（係数ａ、係数ｂ、係数ｃ、係数ｄ）の値が設定されている。

なお流体側補間多項式係数テーブル１２１ｃには、各グリッドにおける流体の速度場の時間微分の計算結果を含めることもできる。
可視化装置１００は、図６と図７に示したシミュレーション結果と、図８に示した事前解析結果Ωとに基づいて、流脈線を計算し可視化する。以下、流脈線の可視化処理について詳細に説明する。

図９は、流脈線可視化処理の手順の一例を示すフローチャートである。以下、図９に示す処理をステップ番号に沿って説明する。
［ステップＳ１０１］流脈線計算部１５０は、事前解析結果記憶部１２０に事前解析ファイル１２１が格納されているか否かを判断する。流脈線計算部１５０は、事前解析ファイル１２１が格納されている場合、処理をステップＳ１０２に進める。また流脈線計算部１５０は、事前解析ファイル１２１が格納されていない場合、処理をステップＳ１０３に進める。

［ステップＳ１０２］流脈線計算部１５０は、事前解析結果記憶部１２０から事前解析ファイル１２１を読み込む。その後、流脈線計算部１５０は、処理をステップＳ１０５に進める。

［ステップＳ１０３］事前解析部１４０は、補間多項式の係数を算出するための事前解析処理を行う。事前解析処理の詳細は後述する（図１０参照）。
［ステップＳ１０４］事前解析部１４０は、事前解析処理によって算出された係数の値を含む事前解析ファイル１２１を生成し、事前解析ファイル１２１を事前解析結果記憶部１２０に格納する。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

次に、事前解析処理（ステップＳ１０３）の詳細について説明する。事前解析処理は、任意時刻ｔにおける心筋の構造情報の推定と血流部分の速度場の推定に使用する係数を生成する処理である。

初めに前提になる出力データについて説明する。出力データには、シミュレーション上の時刻ｔ＝ｔ_iが付与されている。ここで、心筋の構造情報をベクトルＭ（ベクトルｒ，ｔ_i）とする。このとき心筋上には、有限個の離散点ベクトルｒ_kが与えられている。これは、図６に示すように心筋の弾性体情報ファイル１１１ａ，１１１ｂ，・・・に出力されており、離散点はＧＲＩＤ ＩＤで一意に識別され、座標が保存されている。軸ごとの座標の値が、ＧＲＩＤ ＩＤに対応する離散点を示す離散点ベクトルｒ_kの要素である。

有限要素法では、これらの離散点を頂点とした図形（例えば四面体）を要素として、計算が行われる。要素は図６の弾性体情報ファイル１１１ａ，１１１ｂ，・・・のＴＥＴＲＡ ＩＤで一意に識別され、要素を構成する離散点のＧＲＩＤ ＩＤが保存されている。図６の例では四面体要素を用いているが、要素は四面体以外の形状でもよい。これらの離散点の位置と有限要素の配置により心筋の形状が決定される。

流体部分の構造情報も心筋と同様である。流体部分の構造情報をベクトルＢ（ベクトルｒ，ｔ_i）とする。このとき、流体部分も有限個の離散点と有限要素で構成され、図７に示すような流体情報ファイル１１２ａ，１１２ｂ，・・・に出力されている。流体に関して、心筋との違いは、流体情報として、離散グリッド点（ベクトルｒ_k）上での流体の速度場を示す速度場情報（速度ベクトルｖ（ベクトルｒ_k，ｔ_i））も出力されていることである。事前解析部１４０は、すべての時刻ｔ₀，ｔ₁，・・・，ｔ_Lの心筋の構造、血流の構造および速度場情報を用いて、事前解析を実行する。

図１０は、事前解析処理の手順の一例を示すフローチャートである。以下、図１０に示す処理をステップ番号に沿って説明する。
［ステップＳ１２１］事前解析部１４０は、インデックスｉ＝１から順番に、インデックスｉ＝１，２，・・・，Ｎ−１について、ステップＳ１２２のループ処理を繰り返し実行する。

［ステップＳ１２２］事前解析部１４０は、情報読み込み部１３０を介して、時刻ｔ_iのベクトルＢ（ベクトルｒ，ｔ_i）、ベクトルｖ（ベクトルｒ，ｔ_i）、ベクトルＭ（ベクトルｒ，ｔ_i）を、シミュレーション結果記憶部１１０から読み込む。

［ステップＳ１２３］事前解析部１４０は、ステップＳ１２２の処理を１回実行するごとにインデックスｉに１を加算し、インデックスｉ＝Ｎ−１の処理が終了したら、処理をステップＳ１２４に進める。

［ステップＳ１２４］事前解析部１４０は、心筋の構造情報であるベクトルＭ（ベクトルｒ，ｔ_i）を構成する時刻ｔ_iにおけるグリッド座標ベクトルｒ_lの時間微分ベクトルｄｒ_l／ｄｔの値を、すべての時刻、すべてのグリッド点について算出する。

［ステップＳ１２５］事前解析部１４０は、血流の構造情報であるベクトルＢ（ベクトルｒ，ｔ_i）を構成する時刻ｔ_iにおけるグリッド座標ベクトルｒ_lの時間微分ベクトルｄｒ_l／ｄｔの値を、すべての時刻、すべてのグリッド点について算出する。

［ステップＳ１２６］事前解析部１４０は、血流の時刻ｔ_iにおけるグリッド座標ベクトルｒ_l上での速度場ベクトルｖ（ベクトルｒ_l，ｔ_i）の時間微分ベクトルｄｖ（ベクトルｒ_l，ｔ_i）／ｄｔを、すべての時刻、すべてのグリッド点について計算する。

なおステップＳ１２４〜Ｓ１２６に示したインデックスｌは、グリッド番号を示す。また、心筋についてはグリッド座標の総数はＮ_M,elem個存在し、血流についてはグリッド座標の総数はＮ_B,elem個存在するものとする。

次に、心筋、血流部分のグリッド座標ベクトルｒ_lにおける時刻ｔ_iでの時間微分ベクトルｄｒ_l／ｄｔの値を計算方法について詳細に説明する。時間微分ベクトルｄｒ_l／ｄｔの値は、すべてのグリッド点について同じ方法で計算できる。

まず事前解析部１４０は、グリッド番号ｌを指定し、位置ベクトル列ｒ_k（ｔ_i）（ｉ＝０，１，２，・・・，ｎ−１）の情報を取得する。次に事前解析部１４０は、取得した情報からベクトルのＸ，Ｙ，Ｚ成分を抜き出し、座標列Ｘ_l（ｔ_i）（ｉ＝０，１，２，・・・，ｎ−１）、Ｙ_l（ｔ_i）（ｉ＝０，１，２，・・・，ｎ−１）およびＺ_l（ｔ_i）（ｉ＝０，１，２，・・・，ｎ−１）として保存する。ここでｎはファイル数である（ｎは１以上の整数）。

次に事前解析部１４０は、Ｘ成分についてすべての点Ｘ_l（ｔ_i）を通るような滑らかに連続な曲線Ｘ_l（ｔ）を補間法で求める。補間法として、例えばcubic splineを用いる。なお、cubic spline以外の補間法を用いてもよい。cubic splineでは曲線Ｘ_l（ｔ）を区間ｔ_i≦ｔ≦ｔ_(i+1) 毎に定義された補間多項式Ｘ_(l,i)（ｔ）を用いて式（１）で定める。

補間曲線Ｘ_l,i（ｔ）の具体形は

となる。ただし、ｎは心臓の拍動１周期分であるため、事前解析部１４０は、ｎ＋１番目は０番目の値になるように、周期境界条件を課す。例えば事前解析部１４０は、ｔ_n+m＝ｔ_mとする。また事前解析部１４０は、補間曲線Ｘ_l,i（ｔ）の関数値についてもＸ_l,n+m＝Ｓ_l,mとする。ここで係数ａ_l,i，ｂ_l,i，ｃ_l,i，ｄ_l,iは未知数であるが、区間ｔ_i≦ｔ≦ｔ_i+1ごとに微分まで含めて滑らかに連続になる条件を付加することで、連立一次方程式を作ることができる。この連立方程式を数値的に解くことですべての係数ａ_l,i，ｂ_l,i，ｃ_l,i，ｄ_l,iの組を計算することができる（非特許文献５参照）。そして、ｄＸ_l／ｄｔは式（３）から直接計算することができる。

このようにして、時間に関する１階微分および２階微分が連続になるように時刻ｔ_iにおけるｄＸ_l／ｄｔの値を求めることができる。
そして事前解析部１４０は、係数ａ_l,i，ｂ_l,i，ｃ_l,i，ｄ_l,iの配列を事前解析結果Ωに含めて保存する。Ｙ，Ｚ成分についても同様の手続きで補間関数を求めることができる。事前解析部１４０は、このような係数の算出を、すべての心筋のグリッド点と血流のグリッド点について行う。その結果、ベクトルｄｒ_l／ｄｔの値が計算されたこととなる。

次に速度場ベクトルｖ（ベクトルｒ_l，ｔ_i）の時間微分ベクトルｄｖ（ベクトルｒ_l，ｔ_i）／ｄｔを求める手続きを説明する。この時間微分ベクトルｄｖ（ベクトルｒ_l，ｔ_i）／ｄｔは、グリッド点を流体に固定し、グリッド点が流体と共に運動すると仮定した場合におけるグリッド点の加速度を表している。これは流体力学の世界ではラグランジュ座標と呼ばれるものである。

なお、解析対象が非圧縮性流体の場合、加速度に流体の密度を掛けると、そのグリッド点上の粒子に加わる力となる。そのため、グリッド点の加速度を求めることは、非圧縮性流体内のグリッド点上の粒子に加わる力の大きさと方向を求めることに等しい。

流体の速度場は、速度場ベクトルｖ＝（ｖ_x，ｖ_y，ｖ_z）のように成分表示することができる。ここでは速度場のｘ成分（ｖ_x）について説明するが、速度場のｙ成分（ｖ_y）についても同様にして計算することができる。速度場のｚ成分（ｖ_z）を算出する場合、事前解析部１４０は、合成関数の微分を計算する。すると、以下の式（４）が成立する。

式（４）から各グリッド点（グリッド番号ｋ）における時刻ｔ＝ｔ_iでの速度場のｘ成分（ｖ_x）の時間微分を求めることができる。グリッド番号ｋの座標ベクトルｒ_kにおけるベクトルｄｒ_k／ｄｔは、補間多項式からすでに求めている。∇ｖ_xについては、非特許文献６などを参考にして計算できる。残りの速度場の時間に関する偏微分は、以下の式（５）を用いて計算できる。

ここで、ベクトルａ_aveはグリッド番号ｋのグリッド点のラグランジュ座標での平均加速度ベクトル、ベクトルｖ_aveは平均速度ベクトルであり、下記の式（６）、（７）で定義される。

ここで、式（５）にはグリッド点のラグランジュ座標での平均加速度（ベクトルａ_ave）の項と、平均速度（ベクトルｖ_ave）の項が含まれているが、これはグリッド点が任意の運動をしているため、グリッド点自体の運動に起因する見かけの力を補正するためである。

グリッド点が加速度運動をする場合における解析手法として、例えばＡＬＥ（Arbitrary Lagrangian-Eulerian）法がある。
図１１は、ＡＬＥ制御されたグリッド点自身の運動による見かけの力の原理を示す図である。ＡＬＥ制御されたグリッド点は、心臓の収縮に合わせて１拍動すると元の位置に戻り、心臓の拍動に合わせて運動する。このとき、グリッド点は速度が時々刻々と変化するので、加速度が生じる。ところがグリッド点自身に観測装置を設置すれば、観測者は自らを静止したと解釈してしまう。自分自身の加速度運動による効果は見かけの力としてあたかも流体に働くように見えてしまう。そのため、自分自身の加速度運動による補正をしなければならない。この見かけの力を補正するための手続きが式（５）、式（６）、式（７）には含まれている。このようにして、各グリッド点における時刻ｔ＝ｔ_iでの速度場の時間微分を計算することができる。事前解析部１４０は、このような速度場の時間微分をすべてのグリッド点での血流の速度場について求め、各グリッド点での速度場の時間微分を含むデータセットを、事前解析結果Ωとして事前解析結果記憶部１２０に保存する。

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

［ステップＳ１３１］流脈線計算部１５０は、点Ｐ_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を求める。
ステップＳ１４２〜Ｓ１４４における各中間時刻での時間発展は以下のようにして行うことができる。点Ｐ_ijが時刻ｔ＝ｔ_kで、ベクトルｒ_k（ｔ_k）にあったとすると、流脈線の方程式は式（８）で与えられる。

ここでベクトルｖ（ベクトルｒ，ｔ）は時刻ｔ、位置ベクトルｒでの速度（場）である。従って、時間Δｔ後の座標は常微分方程式である式（８）を数値的に解くことによって得られる。式（８）を４次ルンゲ・クッタ法で解くと以下のような計算式が得られる。なお、この考え方はより高次のルンゲ・クッタ型公式にも適用可能であり、高精度の計算が可能である（非特許文献２参照）。ここでは典型例として４次のルンゲ・クッタ公式で説明する。

４次のルンゲ・クッタ公式では、以下の式（９）、式（１０）が成り立つ。

ここで流脈線計算部１５０は、係数α_Iと係数β_Iの値を、ルンゲ・クッタ係数表から取得する。
図１４は、ルンゲ・クッタ係数表の一例を示す図である。図１４に示すルンゲ・クッタ係数表には、Ｉ（Ｉは１〜４の整数）の値ごとの、係数α_Iと係数β_Iの値が示されている。図１４に示すようなルンゲ・クッタ係数表は、例えばメモリ１０２に予め記憶されている。

時刻ｔ＝ｔ_kでの任意の位置ベクトルｒでの速度場ベクトルｖ（ベクトルｒ，ｔ_k）は、ベクトルｖ（ベクトルｒ，ｔ_i）、ベクトルｖ（ベクトルｒ，ｔ_i+1）、ベクトルＢ（ベクトルｒ，ｔ_i）、ベクトルＢ（ベクトルｒ，ｔ_i+1）、ベクトルＭ（ベクトルｒ，ｔ_i）、およびベクトルＭ（ベクトルｒ，ｔ_i+1）から求めることができる。また時刻ｔ＝ｔ_kでの移動境界面をＳ_kとすると、ステップＳ１４２、Ｓ１４３で用いる移動境界面Ｓ_k上の速度場の情報も計算することができる。以下、移動境界面を単に「境界面」と呼ぶ。従って、式（９）、式（１０）、図１４に示したルンゲ・クッタ係数表を用いて時間発展に必要な中間の値を計算することができる。これらの結果を式（８）に代入することで、ベクトルｒ（ｔ_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において、グリッド座標、および速度場は定義されていない。そこで、流脈線計算部１５０は、出力ファイルの速度場から場を近似的に求める。このとき留意すべきは、シミュレーション実行時に、時々刻々とグリッドの位置を動かしていくことである。グリッドの位置の決め方には任意性があるが、心臓など境界自体が移動する問題ではＡＬＥ法がよく利用される。ＡＬＥ法では、シミュレーションに使用する座標自体が、記述する偏微分方程式の解の精度を落とさないように独立に決定される。決定には一般には偏微分方程式が使用される。しかし、データ解析の立場ではグリッドの位置を決定する支配方程式を知ることはできず、与えられたグリッド点の出力値だけが分かっている。このとき、グリッド点の位置が時間に対して連続的に変化する。

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

図１７は、時刻ｔ＝ｔ_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個のグリッド点すべてに対して、ステップＳ１６５〜Ｓ１６８のループ処理を行う。
［ステップＳ１６５］流脈線計算部１５０は、グリッド点ｉが境界面上か否かを判断する。境界面上であれば、流脈線計算部１５０は、処理をステップＳ１６７に進める。境界面上でなければ、流脈線計算部１５０は、処理をステップＳ１６６に進める。

ステップＳ１６５により、グリッド点ｉが境界面上か否かにより、処理が分岐する。このような処理分岐が行われているのは、任意の時刻ｔ_kにおけるグリッド点ｉの座標を補間して求めることになるが、グリッド点ｉが心筋と血流の境界面Ｓ（ｔ）上の点である場合は、血流は心筋に対して滑りなしの拘束条件が課されるためである。

［ステップＳ１６６］流脈線計算部１５０は、グリッド点ｉについての補間式に基づいて、グリッド位置を計算する。補間式は、例えば前述の式（２）である。なお事前解析結果Ωには、グリッド点ｉの座標の時間に関する補間式（式（２））の結果が保存されている。そのため流脈線計算部１５０は、事前解析結果Ωから、任意の時刻におけるグリッド座標を特定することができる。その後、流脈線計算部１５０は、処理をステップＳ１６９に進める。

［ステップＳ１６７］流脈線計算部１５０は、グリッド点ｉが境界面Ｓ（ｔ）上の点であれば、境界条件を満たすように、以下のようにしてグリッド点ｉの位置を計算する。
グリッド点ｉが境界面Ｓ（ｔ）上の点である場合、その速度場から構成される位相空間上の点を（ベクトルｒ_ID（ｔ），ベクトルｖ_ID（ｔ））とする。すると、時刻ｔ_iおよびｔ_i+1では値が出力されているため２点（ベクトルｒ_ID（ｔ_i），ベクトルｖ_ID（ｔ_i））、（ベクトルｒ_ID（ｔ_i+1），ベクトルｖ_ID（ｔ_i+1））を通る曲線として求めればよい。グリッド位置ベクトルｒ_ID（ｔ）と速度場ベクトルｖ_ID（ｔ）について、以下の式（１１）の関係が成立する。

従って、条件式は４つある。そのため、境界面Ｓ（ｔ）上のグリッド点ｉについては、ベクトル_ID（ｔ）を時間について３次式で定めることにより、４つの条件を満たすようにすることができる。従って、ベクトルｒ_ID（ｔ）は、以下の式（１２）から決定できる。

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

となる。ここで、

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

［ステップＳ１６８］流脈線計算部１５０は、補間式（式（１２））を用いて算出した係数を、メモリ１０２に保存する。なお、式（１２）では速度場も同時に求めることができる。速度場についてはここで計算してもよいし、必要に応じて式（１１）および式（１２）を用いて計算することもできる。そのため流脈線計算部１５０は、補間式（式（１２））を用いて算出した係数のみを保存している。

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

場の情報と同様に、心筋の構造情報であるベクトルＭ（ベクトルｒ，ｔ）についても、任意の時刻でのグリッド点情報を計算することができる。
図１８は、時刻ｔ＝ｔ_kにおける弾性体場の情報を求める処理の手順の一例を示すフローチャートである。以下、図１８に示す処理をステップ番号に沿って説明する。

［ステップＳ１７１］流脈線計算部１５０は、読み込まれた心筋の構造情報をメモリに設定する。設定される情報は、ベクトルＢ（ベクトルｒ，ｔ_i）、ベクトルｖ（ベクトルｒ，ｔ_i）、ベクトルＭ（ベクトルｒ，ｔ_i）、ベクトルＢ（ベクトルｒ，ｔ_i+1）、ベクトルｖ（ベクトルｒ，ｔ_i+1）、ベクトルＭ（ベクトルｒ，ｔ_i+1）である。

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

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

［ステップＳ１７５］流脈線計算部１５０は、事前解析結果Ωを用いて補間式（式（２））に従い、グリッド点ｉの位置を計算する。その後、流脈線計算部１５０は、処理をステップＳ１７８に進める。

［ステップＳ１７６］流脈線計算部１５０は、式（１２）に従いグリッド点ｉの座標を計算する。流体におけるすべてのグリッド点ｉの位置を計算することで、流体の構造情報であるベクトルＭ（ベクトルｒ，ｔ_k）を決定することができる。

［ステップＳ１７７］流脈線計算部１５０は、補間式（式（１２））を用いて算出した係数を、メモリ１０２に保存する。
なお流体の処理を先に行った場合、流体側に位置が同じになる対応するグリッド点ｉが存在する。そのため流脈線計算部１５０は、ステップＳ１７６，Ｓ１７７の処理に代えて、ステップＳ１６７における計算結果から、グリッド点ｉの座標を取得し、保存する処理を行ってもよい。

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

次に、任意の時刻ｔにおける血流中のグリッド座標ベクトルｒ_i（ｔ）での速度場ベクトルｖ_l（ベクトルｒ_l（ｔ），ｔ）の計算方法について説明する。
この速度場ベクトルｖ_l（ベクトルｒ_l（ｔ），ｔ）は、流脈線上の点の時間発展をルンゲ・クッタ法で行うための式（９）、式（１０）の計算に用いられる。出力データが与えられている時刻をｔ_iとすると、ｔ_i≦ｔ≦ｔ_i+1を満たす整数ｉを決定することができる。事前解析結果Ωには、各頂点の時刻ｔ_iにおけるグリッド座標ベクトルｒ_l（ｔ_i）の値、その点での速度場ベクトルｖ_l（ベクトルｒ_l（ｔ_i），ｔ_i）の値、および速度場の時間微分ベクトルｄｖ（ベクトルｒ_l（ｔ_i），ｔ_i）／ｄｔが保存されている。従って、速度場は時刻ｔ_iとｔ_i+1で速度場とその微分が与えられていることになる。そのため、微分まで含めて連続な曲線を補間曲線として仮定することができる。与えられた２つの時刻での速度場とその時間微分が与えられているため、式（２０）のような未知変数を４つとする３次の多項式を補間曲線として用いることができる。

係数ベクトルｂ_k,iは、二つの時刻ｔ_iとｔ_i+1において速度場ベクトルｖ_i（ｔ）が滑らかに連続に繋がるように決める。ξ＝ｘ，ｙ，ｚとして、式（２１）〜（２６）を用いて係数を決定することができる。なお、以下の式において、速度の各軸方向の時間微分値を、「ｖ’_x（ｔ）」、「ｖ’_y（ｔ）」、「ｖ’_z（ｔ）」のように「’」を付けて表している。

このようにして、与えられたグリッド点ｉにおける任意の時刻ｔにおける速度場を求めることができる。
次に、任意の時刻ｔにおける血流中の任意の座標ベクトルｒでの速度場ベクトルｖ（ベクトルｒ，ｔ）の計算方法について説明する。

いま、与えられた血流中の任意の座標ベクトルｒは血流中に属するため、座標ベクトルｒを含む血流の有限要素Ｉ_lを求める事ができる。有限要素Ｉ_lに属する頂点の座標をベクトルｒ_l（ｔ）、速度場をベクトルｖ_l（ベクトルｒ_l（ｔ_i），ｔ）とすると、座標ベクトルｒにおける速度場ベクトルｖ（ベクトルｒ，ｔ）は以下の式（２７）から決定することができる。

ここで、Ｎ_i（ベクトルｒ，ｔ）を構造関数と呼ぶ。１次の四面体の構造関数の具体的の計算方法を以下に示す。
図１９は、四面体要素の一例を示す図である。四面体要素Ｉ_lの各頂点の座標が、それぞれＰ₁（ｘ₁，ｙ₁，ｚ₁）、Ｐ₂（ｘ₂，ｙ₂，ｚ₂）、Ｐ₃（ｘ₃，ｙ₃，ｚ₃）、Ｐ₄（ｘ₄，ｙ₄，ｚ₄）で与えられるものとする。このとき、構造関数は式（２８）で与えられる。

ここで、Ｖは四面体要素Ｉ_lの体積である。また、係数ａ_i，ｂ_i，ｃ_iはｉ以外の３つの点が作る平面の方程式の法線ベクトルとして求めることができる。例えば、Ｎ₁は点Ｐ₁以外の３点Ｐ₂，Ｐ₃，Ｐ₄が作る平面の法線ベクトルから決定できる。法線ベクトルをベクトルｎ₁とすると、ベクトルｎ₁＝（ａ₁，ｂ₁，ｃ₁）＝ベクトル（Ｐ₂Ｐ₃）×ベクトル（Ｐ₂Ｐ₄）として計算することができる。ここで、ベクトル（Ｐ₂Ｐ₃）は、点Ｐ₂から点Ｐ₃へのベクトルである。ベクトル（Ｐ₂Ｐ₄）は、点Ｐ₂から点Ｐ₄へのベクトルである。「×」は、ベクトルの外積である。なお、法線ベクトルの向きは底面の３角形Ｐ₂Ｐ₃Ｐ₄からみて、頂点Ｐ₁の向きが正になるようにとるものとする。さらに係数ｄ_iは平面の方程式ａ_iｘ＋ｂ_iｙ＋ｃ_iｚ＋ｄ_i＝０が点Ｐ₂（または点Ｐ₃，点Ｐ₄でもよい）を通ることから決定することができる。従って、任意の点ベクトルｒが四面体要素Ｉ_l中の点であると、式（２８）に従い構造関数の値が決定できる。そして、式（２７）からその場所での速度場ベクトルｖ（ベクトルｒ，ｔ）を計算することができる（非特許文献６参照）。

次に、ルンゲ・クッタ法による時間発展処理について詳細に説明する。
図１３のステップＳ１４２，Ｓ１４３で任意の時刻における心筋と血流のグリッド点ｉの位置を求めた後、ステップＳ１４４で流脈線上の点ベクトルｒ（ｔ_k）の時間発展がルンゲ・クッタ法で行われる。

図２０は、ルンゲ・クッタ法による時間発展処理の手順の一例を示すフローチャートである。以下、図２０に示す処理を、ステップ番号に沿って説明する。
［ステップＳ１８１］流脈線計算部１５０は、インデックスｉ＝１から順番に、インデックスｉ＝１，２，３，４について、ステップＳ１８２〜Ｓ１８３のループ処理を繰り返し実行する。

［ステップＳ１８２］流脈線計算部１５０は、ベクトルｒ（ｔ_k）が属する血流要素番号Ｉ_lを求め、式（２７）による速度ベクトルの計算をするための要素情報を取得する。
［ステップＳ１８３］流脈線計算部１５０は、ルンゲ・クッタ法で用いる中間の速度ベクトルｖ_i（ｉ＝１，２，３，４）を、式（９）、式（１０）で計算する。ここでベクトルｖ_iの計算には、データが与えられていない任意の時刻での流体中の任意の座標における速度場の計算も行われる。例えば、ベクトルｖ₁は時刻ｔ_kにおける座標ベクトルｒ（ｔ_k）における速度場の値であるが、時刻ｔ_kは事前解析によりデータが与えられている時刻とは限らない。また、座標ベクトルｒ_kもデータが与えられている流体中のグリッド点ｉの座標とは限らない。従って、任意の時刻ｔにおける、任意の座標ベクトルｒでの速度場の計算を行うこととなる。

例えば、座標ベクトルｒ（ｔ_k）と時刻ｔ_kが与えられているため、前述した、任意の時刻ｔにおける血流中の任意の座標ベクトルｒでの速度場ベクトルｖ（ベクトルｒ，ｔ）の計算方法（式（２７））により、速度場を計算することができる。なお、ベクトルｖ₂はベクトルｖ₁の情報を含むが、ベクトルｖ₁が既に計算されている。そのため、前述の任意の時刻ｔにおける血流中のグリッド座標ベクトルｒ_i（ｔ）での速度場ベクトルｖ_l（ベクトルｒ_l（ｔ），ｔ）の計算方法（式（２０）および式（２７））でも速度場を計算することができる。

［ステップＳ１８４］流脈線計算部１５０は、ステップＳ１８２〜Ｓ１８３の処理を１回行うごとにｉに１を加算し、ステップＳ１８２〜Ｓ１８３の処理を繰り返す。流脈線計算部１５０は、ｉ＝４に対する計算が終了すると、すべてのベクトルｖ_i（ｉ＝１，２，３，４）の計算が完了したこととなるため、処理をステップＳ１８５に進める。

［ステップＳ１８５］流脈線計算部１５０は、式（９）により、次の時刻ｔ_k+1における流脈線の座標ベクトルｒ（ｔ_k+1）を計算する。
次に、心臓内部での位置判定について詳細に説明する。

ここでは、時間発展の結果得られたベクトルｒ_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］における最大値を数値的に求めることができる。得られた各グリッド点ｉにおける速度場の最大値の中からさらに最大値を選ぶことによって、この二つの時刻の中での速度場のノルム最大の点を最大速度｜ベクトルｖ_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 (7)

1. 流体シミュレーションの解析時間の進行に伴い解析空間内を加速度運動する複数のグリッド点の第１時間間隔の複数の第１時刻の各々における位置を示す位置情報と、前記複数のグリッド点上の流体の前記複数の第１時刻の各々における速度を示す速度情報とを含む流体情報を記憶する記憶部と、
   前記流体情報に基づいて、前記位置情報に表される前記複数のグリッド点の加速度運動に起因する誤差を補正する補正値を含む計算式により、前記複数のグリッド点上の前記流体の前記複数の第１時刻の各々における前記速度の時間微分値をそれぞれ算出し、前記複数のグリッド点上の前記流体の前記複数の第１時刻の各々における前記速度と前記速度の時間微分値とに基づいて、前記解析時間の進行に伴い粒子発生源から順番に出力される複数の粒子の前記第１時間間隔より短い第２時間間隔の複数の第２時刻の各々における位置をそれぞれ算出し、前記複数の粒子の連なりを示す流脈線の表示情報を生成する処理部と、
   を有する流脈線可視化装置。
2. 前記処理部は、前記位置情報に基づいて、前記複数のグリッド点の各々に対し、前記複数の第１時刻の各々におけるグリッド点移動速度とグリッド点移動加速度とをそれぞれ算出し、前記グリッド点移動速度と前記グリッド点移動加速度とを変数として含む前記計算式により、前記複数のグリッド点上の前記流体の前記複数の第１時刻の各々における前記速度の時間微分値をそれぞれ算出する、
   請求項１記載の流脈線可視化装置。
3. 前記処理部は、前記複数のグリッド点の各々について、グリッド点の前記複数の第１時刻の各々における前記流体の前記速度を示す点をそれぞれ接続し、前記点を通るとき、対応する第１時刻での前記流体の前記速度の時間微分値に応じた傾きとなる、前記流体の前記速度の時間変化を表す補間曲線を生成し、前記複数のグリッド点の各々の前記補間曲線に基づいて、前記複数の粒子の前記複数の第２時刻の各々における位置をそれぞれ算出する、
   請求項１または２記載の流脈線可視化装置。
4. 前記記憶部は、さらに、前記解析空間内の構造体の形状の時間変化を示す構造体情報を記憶し、
   前記処理部は、第１解析時刻における第１流脈線に基づいて第２解析時刻における第２流脈線を計算する場合、
   前記第１流脈線上の第１の位置にある離散点を含む、前記解析空間内の部分領域を、前記離散点についての解析対象領域に決定し、
   前記流体情報に示される前記解析対象領域内の前記流体の前記速度に基づいて、前記第２解析時刻における前記離散点上の粒子の移動先を示す第２の位置を計算し、
   前記構造体情報に示される前記解析対象領域内に存在する前記構造体の情報に基づいて、前記第２解析時刻における前記解析対象領域内の前記構造体の占有領域を特定し、
   前記第１の位置と前記第２の位置とに基づいて、前記第２流脈線の前記占有領域への侵入の有無を判定し、
   前記第２流脈線が前記占有領域に侵入しないと判定したとき、前記第２の位置を通過する前記第２流脈線の表示情報を生成する、
   請求項１乃至３のいずれかに記載の流脈線可視化装置。
5. さらに、前記表示情報に基づいて、前記流脈線を表示する表示部を有する、
   請求項１ないし４のいずれかに記載の流脈線可視化装置。
6. コンピュータが、
   流体シミュレーションの解析時間の進行に伴い解析空間内を加速度運動する複数のグリッド点の第１時間間隔の複数の第１時刻の各々における位置を示す位置情報と、前記複数のグリッド点上の流体の前記複数の第１時刻の各々における速度を示す速度情報とを含む流体情報に基づいて、前記位置情報に表される前記複数のグリッド点の加速度運動に起因する誤差を補正する補正値を含む計算式により、前記複数のグリッド点上の前記流体の前記複数の第１時刻の各々における前記速度の時間微分値をそれぞれ算出し、
   前記複数のグリッド点上の前記流体の前記複数の第１時刻の各々における前記速度と前記速度の時間微分値とに基づいて、前記解析時間の進行に伴い粒子発生源から順番に出力される複数の粒子の前記第１時間間隔より短い第２時間間隔の複数の第２時刻の各々における位置をそれぞれ算出し、
   前記複数の粒子の連なりを示す流脈線の表示情報を生成する、
   流脈線可視化方法。
7. コンピュータに、
   流体シミュレーション上の解析時間の進行に伴い解析空間内を加速度運動する複数のグリッド点の第１時間間隔の複数の第１時刻の各々における位置を示す位置情報と、前記複数のグリッド点上の流体の前記複数の第１時刻の各々における速度を示す速度情報とを含む流体情報に基づいて、前記位置情報に表される前記複数のグリッド点の加速度運動に起因する誤差を補正する補正値を含む計算式により、前記複数のグリッド点上の前記流体の前記複数の第１時刻の各々における前記速度の時間微分値をそれぞれ算出し、
   前記複数のグリッド点上の前記流体の前記複数の第１時刻の各々における前記速度と前記速度の時間微分値とに基づいて、前記解析時間の進行に伴い粒子発生源から順番に出力される複数の粒子の前記第１時間間隔より短い第２時間間隔の複数の第２時刻の各々における位置をそれぞれ算出し、
   前記複数の粒子の連なりを示す流脈線の表示情報を生成する、
   処理を実行させる流脈線可視化プログラム。

Images