(57)【要約】【課題】格子生成と形状変更への対応を簡易に行うことができ、かつ精度を保ちつつ計算時間を短縮できる、流体のシミュレーション装置を提供する。【解決手段】流路内において、粘性底層に相当する範囲を特定し、特定された粘性底層に相当する範囲についてプリズム格子を生成し、粘性底層に相当する範囲外では、四面体格子を生成し、当該生成したプリズム格子と四面体格子とを含んだ複合格子を計算用の格子とし、数値流体計算を実行して各格子内での流速、または格子に隣接する流路壁面への剪断応力の少なくとも一方を演算するシミュレーション装置である。【選択図】図１０
2011040055.tif 000002
【特許請求の範囲】
【請求項１】
血管内部の血液の三次元流れをコンピュータに求めさせるための計算格子を作成するプログラムであって、血管壁近くの速度勾配が大きい速度境界層とその影響が強く及ぶ範囲をプリズム押出し格子とし、血管断面の中央部の速度勾配が緩やかな範囲を四面体格子とし、プリズム押出し層の、血管壁に接する第一格子点の距離が摩擦速度と動粘性係数を用いた無次元距離１に相当する距離との差が±２０％の範囲にあることを特徴とするプログラム。
【請求項２】
プリズム押出し層を構成する格子において、血管壁から第二格子点以降の隣接する格子点間距離の比が１．１から１．２５の間にあることを特徴とする請求項１のプログラム。
【請求項３】
プリズム押出し層を形成する血管壁から高さ方向の格子全体の高さが摩擦速度と動粘性係数を用いた無次元距離で１１．６に相当する距離を越える最小の格子数であることを特徴とする請求項２のプログラム。
【請求項４】
血管断面中央部の四面体格子のサイズが、血管断面直径の８％から１０％の間にあることを特徴とする請求項３のプログラム。
【請求項５】
解析の対象となる流路を特定する情報を取得する手段と、
前記取得した情報で特定される流路内において、粘性底層に相当する範囲を演算により特定する手段と、
前記特定された粘性底層に相当する範囲についてプリズム格子を生成する手段と、
前記粘性底層に相当する範囲外では、四面体格子を生成する手段と、
前記生成したプリズム格子と四面体格子とを含んだ複合格子を計算用の格子とし、数値流体計算を実行して各格子内での流速、または格子に隣接する流路壁面への剪断応力の少なくとも一方を演算する手段と、
前記演算された流速または剪断応力の情報を表示する表示手段と、
を含むシミュレーション装置。
【請求項６】
請求項５記載のシミュレーション装置であって、
前記プリズム格子を生成する手段は、流路壁面から最初の格子の位置を表す第１格子位置ｙ0を演算し、当該第１格子位置ｙ0と、流路壁面からの粘性底層高さｙmaxとの間の範囲において、
（ａ）１．１＜ｙi／ｙi-1＜１．２５  （ｉ＝０，１，２…）、
（ｂ）Σｙi＜ｙmax
を満足するｙi（ｉ＝１，２…）を決定し、各格子位置ｙi（ｉ＝０，１，２…）の高さのプリズム格子を生成するシミュレーション装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は血流等を数値流体力学計算によりシミュレートするシミュレーション装置に関する。
【背景技術】
【０００２】
日本国民の死因で循環器系疾患は３０％に及び、今後もさらに増加すると予想されている。患者の血流の情報はせん断応力や圧力など血管負荷を評価することができ、診断・治療において重要である。特許文献１には医用画像から抽出した血管形態をシリコンゴムで再現し、グリセリン水溶液を通流させ血流を予想する例が開示されている。特許文献２には、ボクセル状の簡単な計算格子を用いる一方、超音波ドプラー法による部分的な流速実測結果と計算結果との誤差を仮想的な体積力で仮定して補正し、精度を向上させる技術が開示されている。しかしながら実際には頭部や体内深部にあって超音波ドプラー法が使えないことが多く、信頼性の高い血流シミュレーション方法が望まれている。
【０００３】
数値計算技術とコンピュータ能力の向上により、近年の数値流体力学の進歩は顕著であり、微細な渦や境界層まで直接解析する段階に到達している。その応用範囲は自動車などの産業分野、気象などの自然現象分野のみならず、医療に関わる生体現象分野に及びつつある。血流シミュレーションは直接測定できない血管負荷に関わる物理量を評価することができ、図１のようなこれを用いた診断・治療支援システムの開発は循環系疾患の診断・治療に大きいブレークスルーをもたらすと期待されている。
【０００４】
人体の血管網は種々の口径を持つ多数の血管の分岐と合流を繰り返す複雑な構造であり、数値流体力学解析における離散化に対応する格子分割は極めて重要な役割を持つ。血管流路の３次元容積を格子分割する方法として基本的に２つの方法がある。
【０００５】
一つは六面体格子による構造格子であり、隣接格子間の輸送を担う格子境界面が三次元空間の直行独立した３つの次元によく対応する秩序だった格子分割が可能である。一方、分岐部や合流部の複雑な流路形状に対しては、秩序的な構造を当てはめるのは無理があり歪んだ形の小さい六面体格子を用いる必要がある。このような格子形態の決定には複雑な手順と多大な計算時間を要する。またこのような計算格子では、数値計算の誤差や不安定による発散がしばしば生じることが知られている。
【０００６】
他の一つは四面体格子による非構造格子であり、複雑な形状の流路を均質な四面体格子で容易に分割することができ、また流体計算も比較的速いメリットがある。一方、隣接する格子間の境界面の方向が秩序だっておらず、壁面からの距離も各要素によって異なるため、境界層のように壁面垂直方向に勾配を持つ現象を詳細に解像するには不向きであると言える。
【０００７】
循環系疾患は動脈硬化による血管の閉塞や動脈瘤破裂による出血など血管壁の病変に起因することが多い。血流が血管壁に加える負荷は、圧力とせん断応力であるが、特にせん断応力が血管壁の細胞に大きな影響を及ぼすことが知られている。せん断応力は血管壁に接する速度境界層と呼ばれる速度勾配の大きい薄い層における流動状態が支配的である。
【先行技術文献】
【特許文献】
【０００８】
【特許文献１】特開２００５−４０２９９号公報
【特許文献２】特開２００４−１２１７３５号公報
【発明の概要】
【発明が解決しようとする課題】
【０００９】
血流は基本的に層流であるため、適切な計算分解能があれば特別なモデルを必要とせず六面体格子を用いて速度境界層を精度良く解析することができるが、血管形状のような複雑形状に六面体要素を用いる場合は、計算格子の作成に膨大な時間が必要である。四面体格子は、計算格子の作成は容易であるが、六面体格子と同数の計算格子を用いても、上述した特徴から速度境界層を正確に解析することは困難である。計算分解能を上げた格子数の大きい計算は、計算時間と計算機容量の制約で実行が難しい。
【００１０】
一方、血管壁から遠い血管断面の中心部では速度勾配が緩やかであるため四面体格子を用いても数値拡散がそれほど障害にならない事情がある。計算工学の分野では複雑形状の流れ解析を行う際に壁面近傍をプリズム格子で作成し、その他の部分を四面体格子で作成する手法がしばしば用いられる。しかし、工学の分野における流れ解析と血流解析では、流れの様相が異なるため、工学の分野における流れ解析用の計算格子作成手法をそのまま血流解析用計算格子作成に用いることは難しく、具体的な最適格子分割の指針は確立されていない。
【００１１】
本発明は上記実情に鑑みて為されたもので、格子生成と形状変更への対応を簡易に行うことができ、かつ精度を保ちつつ計算時間を短縮できる、流体のシミュレーション装置を提供することを、その目的の一つとする。
【課題を解決するための手段】
【００１２】
上記従来例の問題点を解決するための本発明は、血管内部の血液の三次元流れをコンピュータに求めさせるための計算格子を作成するプログラムであって、血管壁近くの速度勾配が大きい速度境界層とその影響が強く及ぶ範囲をプリズム押出し格子とし、血管断面の中央部の速度勾配が緩やかな範囲を四面体格子とし、プリズム押出し層の、血管壁に接する第一格子点の距離が摩擦速度と動粘性係数を用いた無次元距離１に相当する距離との差が±２０％の範囲にあることとしたものである。ここで、プリズム押出し層を構成する格子において、血管壁から第二格子点以降の隣接する格子点間距離の比が１．１から１．２５の間にあることとしてもよく、また、プリズム押出し層を形成する血管壁から高さ方向の格子全体の高さが摩擦速度と動粘性係数を用いた無次元距離で１１．６に相当する距離を越える最小の格子数であってもよい。さらに、血管断面中央部の四面体格子のサイズが、血管断面直径の８％から１０％の間にあってもよい。
【００１３】
また、本発明の一態様に係るシミュレーション装置は、解析の対象となる流路を特定する情報を取得する手段と、前記取得した情報で特定される流路内において、粘性底層に相当する範囲を演算により特定する手段と、前記特定された粘性底層に相当する範囲についてプリズム格子を生成する手段と、前記粘性底層に相当する範囲外では、四面体格子を生成する手段と、前記生成したプリズム格子と四面体格子とを含んだ複合格子を計算用の格子とし、数値流体計算を実行して各格子内での流速、または格子に隣接する流路壁面への剪断応力の少なくとも一方を演算する手段と、前記演算された流速または剪断応力の情報を表示する表示手段と、を含む。ここで前記プリズム格子を生成する手段は、流路壁面から最初の格子の位置を表す第１格子位置ｙ0を演算し、当該第１格子位置ｙ0と、流路壁面からの粘性底層高さｙmaxとの間の範囲において、（ａ）１．１＜ｙi／ｙi-1＜１．２５（ｉ＝０，１，２…）、（ｂ）Σｙi＜ｙmaxを満足するｙi（ｉ＝１，２…）を決定し、各格子位置ｙi（ｉ＝０，１，２…）の高さのプリズム格子を生成してもよい。
【発明の効果】
【００１４】
本発明によると、格子生成と形状変更への対応を簡易に行うことができ、かつ精度を保ちつつ計算時間を短縮できる。
【図面の簡単な説明】
【００１５】
【図１】本発明の実施の形態に係るシミュレーション装置の実行する処理の例を表すフローチャート図である。
【図２】本発明の実施の形態に係るシミュレーション装置の処理の対象となる流路の例を表す説明図である。
【図３】本発明の実施の形態に係るシミュレーション装置が生成する格子と、従来の格子との比較例を表す説明図である。
【図４】本発明の実施の形態に係るシミュレーション装置が生成する格子と、従来の格子との比較を表す説明図である。
【図５】本発明の実施の形態に係るシミュレーション装置が生成する格子と、従来の格子とによる解析結果の例を表す説明図である。
【図６】六面体格子分割での格子の例を表す説明図である。
【図７】本発明の実施の形態に係るシミュレーション装置により解析される流路のもう一つの例を表す説明図である。
【図８】本発明の実施の形態に係るシミュレーション装置により生成される複合格子の例を表す説明図である。
【図９】本発明の実施の形態に係るシミュレーション装置が生成する格子と、従来の格子とによる解析結果の例を表す説明図である。
【図１０】本発明の実施の形態に係るシミュレーション装置の構成例を表すブロック図である。
【図１１】本発明の実施の形態に係るシミュレーション装置により解析される流路のさらにもう一つの例を表す説明図である。
【図１２】本発明の実施の形態に係るシミュレーション装置が生成する格子と、従来の格子とによる解析結果の例を表す説明図である。
【発明を実施するための形態】
【００１６】
［概要］
血管形状の計算格子作成にあたっては、壁面せん断応力（壁面近傍の速度境界層）を精度よく解像可能であり、かつ管内の流速分布を精度よく解像可能である必要がある。また、その上で計算負荷を抑制する為に、計算格子数は最小限に留める必要がある。血流の解析においては、管内の流れは、ほとんどの場合層流である事や、解析可能な血管の直径などがある一定の範囲内であるため、計算格子を作成するにあたって、ある一定の基準が作成可能であると考えられる。
【００１７】
後に述べる本実施の形態において解決されるべき課題の一つは、せん断応力を支配する速度境界層の計算精度の確保と計算処理速度向上の両立を達成可能とするために、プリズム格子と四面体格子をそれぞれどの領域にどのような寸法でどの程度配置すればよいかの基準を、血流解析以前の既知の情報を基に決定する事である。
その際に、（１）プリズム層の第１格子点距離、（２）プリズム層の拡張比、（３）プリズム層の総合高さ、（４）四面体格子一辺の長さの基準が必要となる。
【００１８】
速度境界層を十分な精度で解像可能な計算格子を作成するにあたっては、以下のように表わされる無次元距離y+が１程度となるように壁面からの第１格子点（壁面に最も近い格子点）をとる事が一つの基準とされている。
【数１】
2011040055.tif 000003
ここで
【数２】
2011040055.tif 000004
であり、yminは壁面からの第１格子点の距離（第１格子点距離）、u*は壁面摩擦速度、νは動粘性係数である。これら（１）、（２）式と血管内の流れが層流であるという仮定から、血流解析用の計算格子の第１格子点距離の基準を以下のように導く。
【００１９】
すなわち、ハーゲン・ポアズイユ流れの流速の式から血管半径と壁面摩擦速度の関係を導くと、第１格子点の距離yminは以下のようになる。
【数３】
2011040055.tif 000005
ここで
【数４】
2011040055.tif 000006
である。y+=1とおき、血液の平均的な密度、粘性、及び最大流速を（３），（４）式に与えると、第１格子点の距離yminが見積もられる。見積もられた値を中心にプラスマイナス２０％程度の範囲を最適範囲とする。この最適範囲を超えてyminを大きく取り過ぎると壁面近傍での境界層を十分な精度で解析出来なくなる。また、最適範囲を下回ってyminを小さく取り過ぎると、プリズム要素を構成する３角形の辺とプリズム格子高さのアスペクト比が大きくなり過ぎるため、壁面近傍の計算格子の品質が低下する。
【００２０】
壁面垂直方向の隣接格子間隔比の基準値としては、一般的に1.2程度以内に収めるのが望ましいとされる。血流解析用の計算格子作成にあたっては、拡張比の範囲を1.1以上、1.25以下とし、第１格子点の距離と後述するプリズム層の総合高さを満たすようにこの範囲で拡張比を調整してプリズム層を作成する。
【００２１】
境界層を十分な精度で解像するために、境界層内の粘性底層と呼ばれる壁面近傍の速度勾配が大きい領域内にプリズム層を配置出来るようにプリズム層の総合高さを決定する。粘性底層はy+<11.6のように定義されるのが一般的である。（３）式でy+=11.6とすれば粘性底層の厚さを求めることが出来る。動脈瘤の好発部位であるWillis動脈輪のように、解析の対象となる血管形状は、血管直径が２mm〜２０mm程度であり、流入・流出部の口径比が３以下である場合が多い。これらの事に併せて中心部分の四面体格子との接続等も考慮して試作を行った結果、プリズム層総合高さを血管形状の最小断面直径の２５％程度（±５％）とした時に、全体の品質が良く、血管全体にわたって粘性底層を効率良く解像できる計算格子が作成可能である。
【００２２】
四面体格子一辺の長さについては、上述したプリズム層との接続や、全体のバランスを考慮して試作を行った結果、血管形状の流入・流出断面においてはその直径の１０％程度とし、血管壁面上では最大断面直径の８％〜１０％程度とした時に全体のバランスが良い計算格子が作成可能であった。四面体格子一辺の長さを小さく取り過ぎると、計算負荷が膨大になるにも関わらず、得られる計算結果は最適値の結果とほとんど同じものとなってしまう。大きく取り過ぎると、血管内の流れを十分な精度で解像することが出来ない。
【００２３】
このようにすると、血流シミュレーションにおいて血管断面の流れの状況および数値拡散の障害の程度に応じて格子の種類と配置を配分し、格子形成および計算実行に要する時間を節減することができる。また数値的に計算格子の種類と配置と寸法を指定するので、医用画像から計算格子生成プログラムの自動化も容易である。
【００２４】
［ハードウェア］
本発明の実施の形態に係るシミュレーション装置について、図面を参照しつつ説明する。なお、本実施の形態のシミュレーション装置は、血管内の血液の流れ（血流）をシミュレートする場合だけでなく、三次元空間内に規定された流路を流れる流体全般に適用可能なものであるが、以下の例では血流をシミュレートする場合を例として説明する。
【００２５】
図１０は、本実施の形態に係るシミュレーション装置１の一例を表す構成ブロック図である。図１０に例示するように、本実施の形態のシミュレーション装置１は、制御部１１と、記憶部１２と、操作部１３と、表示部１４とを含んで構成されている。
ここで制御部１１は、ＣＰＵなどのプログラム制御デバイスであり、記憶部１２に格納されたプログラムに従って動作する。本実施の形態の制御部１１は、まずシミュレーションの対象となる流体の流路を特定する情報を取得する。そして当該流路内に、数値流体計算用の格子を生成する。制御部１１は、この格子を用いて数値流体計算を実行し、その計算の結果を出力する。本実施の形態において特徴的なことの一つは、この制御部１１が生成する格子が、プリズム格子と、四面体格子とを組み合わせたものとなっていることである。ここでプリズム格子は、略三角柱状をなす格子であり、流路の形状によっては、この格子を隙間なく配するために、対向する三角形の面の一方が、他方に比べて小さくなるようテーパーをつける。この制御部１１の詳しい処理の内容については後述する。
【００２６】
記憶部１２は、メモリデバイスや、ディスクデバイスなどにより実現される。この記憶部１２は、制御部１１によって実行されるプログラムを保持している。本実施の形態では、このプログラムは、ＤＶＤ−ＲＯＭなどのコンピュータ可読な記録媒体に格納されて提供され、この記憶部１２にインストールされたものであってもよい。また、記憶部１２は、制御部１１のワークメモリとしても動作する。
【００２７】
操作部１３は、キーボードやマウス等であり、利用者の指示操作を受け入れて制御部１１に出力する。表示部１４は、ディスプレイ等であり、制御部１１から入力される指示に従って情報を表示出力する。また、このシミュレーション装置１は、ネットワークインタフェースなどを備えてもよい。ネットワークインタフェースを備える場合、シミュレーション装置１は、上記流体の流路を特定する情報を、ネットワークを介して取得することとしてもよい。
【００２８】
［制御部１１の処理］
ここで、制御部１１による具体的な処理について述べる。この制御部１１は、図１に例示するように、ＣＴ画像（Computed Tomography image）や、ＭＲＩ（Magnetic Resonance Imager）の画像の入力を受け入れる（Ｓ１）。そして制御部１１は、これらの画像を分析して三次元的な流路の情報（血管形状の情報）を抽出する（Ｓ２）。このようにＣＴ画像や、ＭＲＩの画像に基づいて三次元的な血管形状の情報を得る方法は、広く知られているので、ここでの詳細な説明を省略する。なお、ここで三次元的な流路の情報は、例えば流路をメッシュで表現し、当該メッシュの格子点の三次元座標として表すこととしておけばよい。
【００２９】
制御部１１は、抽出した流路の情報を参照し、流路内に数値流体計算用の格子を生成する（Ｓ３）。流路の内壁面（血管壁、以下壁面と呼ぶ）のごく近傍では、分子粘性の作用が支配的になる粘性底層が生じ、粘性底層の外側（壁面から粘性底層の厚さより離れたところ、血管の中心側）では流速が滑面対数分布則に従い、速度勾配が緩やかであるような層ができている。
【００３０】
本実施の形態の制御部１１は、この粘性底層には三角柱状のプリズム格子を隙間なく配し、また、粘性底層よりも血管の中心側の範囲については、四面体格子を隙間なく配した計算用の格子を生成する。この格子の概要を図示すると、図８に示すようになる。すなわち、管面を平面視したときには、三角形が配列された形状となり（図８（ａ））、管面を流路に垂直な面で切った断面でみると、粘性底層には略三角柱状を側面から見たものとなり、四角形状が配列された状態となる。また、粘性底層より内側（血管中心側）では、四面体形状が配列された状態を見ることになる（図８（ｂ））。
【００３１】
なお、これらの各層では、流速を表す式が異なるのであるが、その境界部では両者の速度は一致することが当然に求められる。そこで、これらの流速の式が互いに等しいとして方程式を生成し、これを（１）式の無次元距離y+について解くと、y+＝11.6という値が得られる。すなわち、粘性底層の厚さは、y+＝11.6として表すことができる。
【００３２】
制御部１１による上記格子の生成方法は広く知られた方法を採用できるが、その一例を次に述べる。制御部１１は、複数のプリズム格子２１の一方の底面（テーパーをかける場合は大きい方の三角形の面）を壁面２０に接するようにして隙間なく配する。このように、任意の面上に三角形を隙間なく配する方法については広く知られた方法を採用できる。
【００３３】
ここで制御部１１は、第１格子点距離ｙminを、次のようにして計算する。以下では格子内の局所的な血流をハーゲン・ポアズイユ流れ（Hagen-Poiseuile flow）であると仮定しておく。壁面２０に配した三角形の頂点の一つを注目点として、この注目点を含み、血管の軸（流れの方向）に垂直な面において、血管の内面（壁面）の平均半径をａとする。また、この面における血管断面を円で近似したときの円の中心（この円の半径が平均半径となる）からの距離をｒととり、壁面摩擦速度u*を
【００３４】
【数２】
2011040055.tif 000007
で定義する。
【００３５】
なお、νは対象となる流体（ここでは血液）の動粘性係数である。そして制御部１１は、壁面からの第１格子点の距離を、
【００３６】
【数１】
2011040055.tif 000008
の無次元距離y+を１とするようなyminとして定める。
【００３７】
具体的に、ハーゲン・ポアズイユ流れの流速の式を用い、血管半径と壁面摩擦速度の関係を導くと、第１格子点の距離yminは（１）式より、以下のようになる。
【００３８】
【数３】
2011040055.tif 000009
ここで
【数４】
2011040055.tif 000010
である。制御部１１は、この（３），（４）式を用い、y+=1とおいて、血液の平均的な密度、粘性、及び最大流速（既知である）を代入して、第１格子点の距離yminを得る。
【００３９】
なお、ここでは制御部１１は、このようにして算出した第１格子点の距離yminを用いて、以下の処理を行うものとするが、実験的に、この値を中心にプラスマイナス２０％程度の範囲であればよいことが知られている。
【００４０】
すなわち、制御部１１は、他の要因により、ここで算出した第１格子点の距離yminを仮の距離として、この仮の距離に基づき、実際に計算に用いる第１格子点の距離y′minを算出してもよい。この際、y′min／yminが０．８から１．２の間にあるようにすればよい。
【００４１】
なお、既に述べたように、上記最適範囲を超えてyminを大きく取り過ぎると壁面近傍での境界層を十分な精度で解析出来なくなる。また、最適範囲を下回ってyminを小さく取り過ぎると、プリズム要素を構成する３角形の辺とプリズム格子高さのアスペクト比が大きくなり過ぎるため、壁面近傍の計算格子の品質が低下する。
【００４２】
次に制御部１１は、プリズム格子を配する範囲（壁面からの距離）を決定する。既に述べたように、粘性底層はy+<11.6の範囲にあるから、制御部１１は、（３）に相当する
【００４３】
【数５】
2011040055.tif 000011
を用いてy+=11.6を代入し、粘性底層の厚さymaxを求める。
【００４４】
制御部１１は、予め定められた隣接格子間隔比（ここでの例のように血液を対象とする場合は経験的に1.1以上、1.25以下とすればよい）に基づいて、第１格子点の距離y0=yminとし、また以下の格子点（第２，３…格子点）までの距離y1，y2，…を、
（ａ）１．１＜ｙi／ｙi-1＜１．２５  （ｉ＝０，１，２…）、
（ｂ）Σｙi＜ｙmax
の各条件を満足するように定める。ここで（ｂ）ではｉについて総和する。
【００４５】
例えば、第１格子点の距離y0=yminとし、以下、上記（ａ）の条件を満足するよう、順次y1，y2，…を、（ｂ）の条件を満足する間だけ繰り返して求めていけばよい。
【００４６】
なお、ｙi／ｙi-1の比は、Σｙiが最もｙmaxに近くなる値として求めればよい。制御部１１のようなＣＰＵなどにより、このような値を求める方法は広く知られている方法があるので、ここでの詳しい説明を省略する。
【００４７】
なお、粘性底層の厚さは血管形状や、血管の径、流入・流出部の口径比等によっても影響される。そこで、プリズム格子を配する範囲（壁面からの距離）を予め、処理Ｓ２にて求めた血管形状の中で断面直径が最小となる径ｒminを見出し、当該のｒminの２５％（±５％）に相当する無次元距離y+を、ｙmaxとしてもよい。
【００４８】
制御部１１は、壁面から、以上のように求めた第ｉ格子点の距離ｙi（ｉ＝０，１，２…）の位置にそれぞれ仮想の壁面をつくっておく。制御部１１は、壁面を互いに隣接した三角形からなるメッシュに分割する（三角分割）。そしてこの三角分割で得られる各三角形の頂点座標Ｔ1，Ｔ2…を得ておく。
【００４９】
次に制御部１１は、第１格子点の距離の位置に生成した仮想の壁面においても同じ方法で三角分割を行い、壁面を三角分割して得られた各三角形の頂点座標Ｔ1，Ｔ2…に対応する各三角形の頂点座標Ｔ0_1，Ｔ0_2…を得る。同様に、第２、第３…格子点（第ｉ格子点）の距離の位置に生成した仮想の壁面においても同じように、壁面を三角分割して得られた各三角形の頂点座標Ｔ1，Ｔ2…に対応する各三角形の頂点座標Ｔi_1，Ｔi_2…を得る。
【００５０】
そして、壁面及び各仮想的な壁面において対応する頂点座標同士を結び、稠密に配された、略三角柱の格子からなるプリズム層を形成する。
【００５１】
制御部１１は、最も血管の中心側にある仮想的な壁面において、三角分割により形成した三角形を面の一つとした四面体を含んだ四面体格子２２を、血管の中心側、プリズム層を形成していない範囲に稠密に配する。この方法は、広く知られた四面体格子への分割処理を採用できるので、ここでの詳しい説明は省略する。
【００５２】
なお、制御部１１は、血管内面である壁面２０における三角分割の際に、当該壁面２０上に生成する三角形の一辺のサイズを次のようにして定めてもよい。すなわち、四面体格子２２の一辺の長さについては、上述したプリズム層との接続や、全体のバランスを考慮しつつ行った実験により、血管形状の流入・流出断面においてはその直径の１０％程度、また血管壁面上では最大断面直径の８％〜１０％程度とすることが適当であることが分かった。これは、四面体格子２２の一辺の長さを上記の程度より小さく取り過ぎると、計算負荷が膨大になるにも関わらず計算精度の向上が期待できず、また上記の程度より大きく取り過ぎると、血管内の流れを十分な精度で解像できないことに基づく。
【００５３】
従って、これにより、例えば四面体格子２２の一辺の長さを血管断面の１０％としておくと、最も血管の中心側にある仮想的な壁面において、三角分割により形成した三角形の一辺の長さもまた、血管断面の１０％となるから、当該最も血管の中心側にある仮想的な壁面の中心からの距離と、壁面（血管の内面）までの中心からの距離の比とに基づいて、壁面（血管の内面）２０上に生成する三角形の一辺のサイズを定めることができる。
【００５４】
なお、既に説明したように、制御部１１は、他の方法により、粘性底層には三角柱状のプリズム格子２１を隙間なく配し、また、粘性底層よりも血管の中心側の範囲については、四面体格子２２を隙間なく配した計算用の格子を生成することとしてもよい。
【００５５】
制御部１１は、以上の処理により、数値流体計算用の計算格子を生成すると、この計算格子を利用して流体のシミュレーション処理を実行する（Ｓ４）。このシミュレーション処理については、広く知られた数値中退計算方法を用いることができるので、詳しい説明を省略するが、このシミュレーション処理により、各計算格子内の流体の流速の時間変化が計算される。
【００５６】
制御部１１は、各計算格子内の流体流速の時間変化を得たところで、各計算格子に隣接する壁面（血管内面）に対する応力を演算する（Ｓ５）。そして制御部１１は、表示部１４に対して処理Ｓ２にて得た三次元の流路の情報に基づいて生成した三次元グラフィクス（このグラフィクスの生成もまた、広く知られているのでその詳細な説明は省略する）において、上記計算格子に隣接する壁面ごとに演算された当該応力の値に関連づけて予め定められた色の情報を取得し、上記三次元グラフィクスにおいて、対応する計算格子に隣接する壁面の部分を、当該取得した色に着色する（Ｓ６）。これにより得られた画像の例を図９（ａ）に示す。
【００５７】
利用者は、操作部１３から指示を入力し、この画像を三次元的に回転等させつつ、応力の集中する位置などを確認するなどの作業を行うことになる。
【実施例】
【００５８】
［実施例１］
基本的な血管形態である分岐部を模擬したものを図２に示す。この図２に示したＴ字分岐管を、六面体格子に分割した場合と、四面体格子に分割した場合と、本実施の形態の制御部１１により、プリズム格子と四面体格子との組み合わせ（以下、複合格子と呼ぶ）に分割した場合とのそれぞれの処理時間（計算格子の作成時間）を図３に示す。なお、ここで解像度は略同じとした。この結果、本実施の形態による複合格子の生成時間は六面体格子の数百分の１であり、四面体格子とほぼ同じであった。
【００５９】
さらに、これらの計算格子を用いて数値流体計算を行ったときの、計算時間の比較を図４に示す。本実施の形態の複合格子では、後に述べるように、タイムステップを大きく（ステップ数を小さく）できるため、六面体格子より格子数を増加して解像度を上げたにもかかわらず、計算時間は、六面体格子の１／１０以下に節減できた。なお図４では複合格子においてタイムステップを大きく（ステップ数を小さく）しない例も、比較のために含めている。
【００６０】
せん断応力評価のベースとなる壁面付近の速度分布の計算結果を図５に比較する。本実施の形態の複合格子は六面体格子を用いた高精度計算とよく一致する境界層速度分布を得ることができた。なお、六面体格子分割では、図６のように分岐管接続部で歪の大きい小さい格子が生成され、タイムステップに制約が生じやすい。
【００６１】
具体的に六面体格子のクーラン数Ｃ（流速ｕ、時間刻み幅Δｔ、格子幅Δｘの関係式で、Ｃ＝ｕΔｔ／Δｘ）は、２．２９８であるのに対し、本実施の形態による複合格子のクーラン数Ｃは、１．４４２であり、流速が同じ場合であれば、Δｔ／Δｘを六面体格子より小さくできる。
【００６２】
［実施例２］
さらに別の実施例として、図７に示す医用画像から得た頚動脈形態を用いた血流シミュレーションの解析条件を図８に示す。本実施の形態の複合格子は従来の四面体格子と比較して図９のように滑らかで妥当なせん断応力分布を計算することができた。四面体格子は秩序だった格子は位置ではないため、分布計算にむらが生じやすい。
【００６３】
なお、図８の計算では、離散化手法として、
空間：有限体積法（対流項：2次精度中心差分）
時間：完全陰解法、時間刻み幅Δ t = 1.0×10-3秒
とし、境界条件として、
流入：0.4m/s (一様)、
壁面：スリップなし、剛体壁、
流出：自由流出境界条件、
とした。また計算格子（解析格子）は、本実施の形態の複合格子とした。このとき要素数は225,114であった。
【００６４】
六面体格子や本実施の形態の複合格子においては、格子が壁面垂直方向に直交するものであるため、壁面にかかるせん断応力の連続性が維持されやすい。
【００６５】
［実施例３］
次に、流路形状を中大脳動脈（ＭＣＡ：図１１参照）とし、流入条件を平均流速０．２８ｍ／ｓのポアズイユ流れとし、流出条件を自由流出であるとし、レイノルズ数を５７５とおいたときの数値流体計算を行った結果について説明する。
【００６６】
ここで数値流体計算の計算スキームは２次中心差分、時間積分にはクランク−ニコルソン法を採用した。また、結果の出力にあたっては、０．０２秒間隔で平均化したものを用いた。
【００６７】
本実施形態の複合格子、六面体格子、四面体格子のそれぞれに対して、上記数値流体計算を行った結果は、図１２（ａ）から（ｃ）に示すようなものとなった。図１２（ａ）は、本実施形態の複合格子を利用した結果であり、図１２（ｂ）は六面体格子を利用した結果である。これら本実施形態の複合格子を利用した結果と、六面体格子を利用した結果（せん断応力の分布の表示結果）とは、ほとんど同じものとなった。一方、図１２（ｃ）に示した四面体格子を利用した結果では、せん断応力の分布に連続性がなく、むらが生じたり、応力が集中するべき位置において応力集中の状態が明瞭に現れないなどの現象が見られた。
【００６８】
さらに、計算ステップ（時間刻み幅Δｔ）を、１．０ｄ−４、５．０ｄ−４、１．０ｄ−３と変化させて計算させたときの計算時間（時間：分）を次の表に示す。
【００６９】
【表１】
2011040055.tif 000012
【００７０】
なお、六面体格子は、時間刻み幅Δｔ＝１．０ｄ−３のとき、計算結果が収束せずに発散した。また、四面体格子のうちΔｘが比較的粗いもの（coarse）については、流れが再現できず、妥当な結果が得られなかった。
【００７１】
この表に示すように、本実施の形態の複合格子では、時間刻み幅Δｔ＝１．０ｄ−３においても発散せず、２８時間１８分にて計算が終了し、他の格子に比べてより短時間で妥当な結果が得られる。
【００７２】
また計算格子を生成する時間を比較すると、次の表のようになった。
【００７３】
【表２】
2011040055.tif 000013
【００７４】
以上のように、本実施の形態のシミュレーション装置では、粘性底層に相当する部分を略三角柱状のプリズム格子とし、他の速度勾配が緩やかであるような層に相当する部分を四面体格子とした複合格子を生成して、数値流体計算に用いるので、格子作成時間や計算時間は比較的四面体格子に近いため格子生成と形状変更への対応を簡易に行うことができ、六面体格子に比べて低減された時間となり、解析の結果は六面体格子に近く、連続した妥当な結果が得られる。すなわち、精度を保ちつつ計算時間を短縮できる。
【００７５】
なお、本実施の形態のシミュレーション装置は、循環系疾患の診断、治療、医療機器の設計等への応用が可能であるばかりでなく、一般の流体計算にも応用が可能である。
【符号の説明】
【００７６】
１  シミュレーション装置、１１  制御部、１２  記憶部、１３  操作部、１４  表示部、２０  壁面、２１  プリズム格子、２２  四面体格子。
【図１】
2011040055.tif 000014
【図２】
2011040055.tif 000015
【図３】
2011040055.tif 000016
【図４】
2011040055.tif 000017
【図５】
2011040055.tif 000018
【図６】
2011040055.tif 000019
【図７】
2011040055.tif 000020
【図８】
2011040055.tif 000021
【図９】
2011040055.tif 000022
【図１０】
2011040055.tif 000023
【図１１】
2011040055.tif 000024
【図１２】
2011040055.tif 000025