(57)【要約】      （修正有）【課題】動脈瘤が成長するか否かを精度よく判定することができる方法を提供する。【解決手段】動脈瘤に流入する血液の流入エネルギーＥに基づいて動脈瘤が成長するか否かを判定する。この方法は、動脈瘤に流入する血液の密度をρ、動脈瘤に接続している血管と動脈瘤との境界平面に直交する方向における動脈瘤に流入する血液の流速を｜ｖ’in｜、動脈瘤に流入する血液の体積流量をＱin、動脈瘤に流入する血液の圧力をＰin、とするとき、流入エネルギーＥを、Ｅ＝（１／２）・ρ・Ｑin・｜ｖ’in｜２＋Ｐin・Ｑinなる式、又は、Ｅ＝Ｐin・Ｑinなる式、に基づいて算出する。【選択図】図１０
2015073858.tif 000002
【特許請求の範囲】
【請求項１】
動脈瘤が成長するか否かを判定するための動脈瘤成長判定方法であって、
前記動脈瘤に流入する血液の流速及び圧力に基づいて同動脈瘤に流入する同血液の流入エネルギーＥを求め、
前記流入エネルギーＥに基づいて前記動脈瘤が成長するか否かを判定する、
方法。
【請求項２】
請求項１に記載の動脈瘤成長判定方法において、
前記流入エネルギーＥを求めるステップは、
前記流入エネルギーＥを、前記動脈瘤に流入する血液の流速に基づいて得られる体積流量Ｑinと、前記動脈瘤に流入する血液の圧力Ｐinと、の積Ｐin・Ｑinである圧力エネルギーに基づいて求めるステップである、
方法。
【請求項３】
請求項２に記載の動脈瘤成長判定方法において、
前記流入エネルギーＥを求めるステップは、
前記動脈瘤に流入する血液の密度をρ、
前記動脈瘤に接続している血管と同動脈瘤との境界平面に直交する方向における前記動脈瘤に流入する血液の流速の大きさを｜ｖ’in｜、
とするとき、
前記流入エネルギーＥを、
Ｅ＝（１／２）・ρ・Ｑin・｜ｖ’in｜２＋Ｐin・Ｑin
なる式に基づいて求めるステップである、
方法。
【請求項４】
請求項１乃至請求項３の何れか一項に記載の動脈瘤成長判定方法に使用される装置であって、
前記動脈瘤及び同動脈瘤に接続した血管を非圧縮性流体としての血液が流れる構造体として扱い、前記血管の内部を流れる前記血液に有限体積法を適用することにより前記血液の流速及び圧力を流体解析により取得し、前記取得した流速及び圧力を用いて前記流入エネルギーＥを求める装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、動脈瘤が成長するか否かを判定するための方法及びその方法に使用される装置に関する。
【背景技術】
【０００２】
従来から知られている動脈瘤の有無を判定する方法の一つは、動脈瘤の発生が予想される部分を経た血管経路の脈波信号と動脈瘤の発生が予想される部分とは異なる部分を経た血管経路の脈波信号とを検出し、それらの脈波信号の周波数特性を比較することによって動脈瘤の有無を判定する（例えば、特許文献１を参照。）。この他にも、動脈瘤の有無を判定する方法が幾つか提案されている（例えば、特許文献２及び３を参照。）。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特開２０１３−９４２６４号公報
【特許文献２】国際公開第ＷＯ２００９／０５０３０４号
【特許文献３】国際公開第ＷＯ００／０５３０８１号
【発明の概要】
【０００４】
上記従来の技術は動脈瘤の有無を判定する方法である。しかしながら、動脈瘤が成長がするか否かを判定する方法は殆ど知られていない。更に、動脈瘤の成長には種々の因子が関連しているので、動脈瘤が成長するか否かを正確に判定することは困難である。そのため、動脈瘤が成長するか否かを精度よく判定することができる方法が求められている。
【０００５】
本発明は動脈瘤が成長するか否かを精度よく判定することができる方法を提供する。
より具体的に述べると、本発明の動脈瘤成長判定方法は、
第１ステップとして、前記動脈瘤に流入する血液の流速及び圧力に基づいて同動脈瘤に流入する同血液の流入エネルギーＥを求め、
第２ステップとして、前記流入エネルギーＥに基づいて前記動脈瘤が成長するか否かを判定する、
ステップを含む方法である。
【０００６】
この場合、前記流入エネルギーＥを求めるステップは、
前記流入エネルギーＥを、「前記動脈瘤に流入する血液の流速に基づいて得られる体積流量Ｑinと、前記動脈瘤に流入する血液の圧力Ｐinと、の積Ｐin・Ｑin」である圧力エネルギーに基づいて求めるステップであってもよい。
【０００７】
更に、前記流入エネルギーＥは、例えば、
Ｅ＝（１／２）・ρ・Ｑin・｜ｖ’in｜２＋Ｐin・Ｑin
なる式に基づいて求めることができる。
ここで、
ρは前記動脈瘤に流入する血液の密度、
｜ｖ’in｜は前記動脈瘤に接続している血管と同動脈瘤との境界平面に直交する方向における前記動脈瘤に流入する血液の流速（の大きさ）、
Ｑinは前記動脈瘤に流入する血液の体積流量、
Ｐinは前記動脈瘤に流入する血液の圧力、
である。
【０００８】
即ち、流入エネルギーＥは、前記境界平面を通って動脈瘤に流入する血液の「運動エネルギー（（１／２）・ρ・Ｑin・｜ｖ’in｜２）と圧力エネルギー（Ｐin・Ｑin）との和」として求めることができる。但し、流入エネルギーＥは、前記圧力エネルギー（Ｐin・Ｑin）のみを含むエネルギーであってもよい。
【０００９】
発明者は、複数の患者のそれぞれについての前記流入エネルギーＥと動脈瘤の表面積増加量との関係を調査した。その結果を図９及び図１０に示す。図９は、前記流入エネルギーＥとして前述した血液の「圧力エネルギーと運動エネルギーとの和」を採用したグラフである。図１０は、前記流入エネルギーＥとして前述した血液の「圧力エネルギー」を採用したグラフである。表面積増加量とは、各患者の一年あたりの動脈瘤の表面積の増加分である。
【００１０】
図９及び図１０から、流入エネルギーＥが大きいほど動脈瘤の表面積増加量が大きいことが理解される。即ち、流入エネルギーＥが大きいほど動脈瘤が成長することが判明した。従って、流入エネルギーＥに基づいて動脈瘤が成長するか否かを精度良く判定することができる。
【００１１】
より具体的な一つの例においては、前記流入エネルギーＥが所定の閾値以上であるか否かを判定し、前記流入エネルギーＥが前記閾値以上であると判定された場合に動脈瘤が成長すると判定する。これによれば、動脈瘤が成長するか否かを簡単に且つ精度よく判定することができる。勿論、本発明は係る具体例に限定されることはなく、流入エネルギーＥが動脈瘤の成長に相関が強いパラメータであることを用いた動脈瘤成長判定方法であれば、どのような方法であっても本発明に含まれる。
【００１２】
更に、本発明の他の一つは、前記動脈瘤成長判定方法に使用される装置である。この装置は、前記動脈瘤及び同動脈瘤に接続した血管を「非圧縮性流体としての血液」が流れる構造体として扱い、前記血管の内部を流れる前記血液に有限体積法を適用することにより前記血液の流速及び圧力を流体解析により取得し、前記取得した流速及び圧力を用いて前記流入エネルギーＥを求める装置である。この装置によれば、前記流入エネルギーＥを容易に取得することができる。
【００１３】
本発明の他の目的、他の特徴及び付随する利点は、以下の図面を参照しつつ記述される本発明の実施形態についての説明から容易に理解されるであろう。
【図面の簡単な説明】
【００１４】
【図１】ＣＴスキャン（又はＭＲＩ）を用いて撮像された「動脈瘤及びその周辺部」の３次元画像を模式的に示した図である。
【図２】流入エネルギーを算出するための境界条件を説明するための図である。
【図３】流入エネルギーの算出方法を説明するための図である。
【図４】ある患者の動脈瘤の表面積及び流入エネルギーの経年変化を示すグラフである。
【図５】別の患者の動脈瘤の表面積及び流入エネルギーの経年変化を示すグラフである。
【図６】別の患者の動脈瘤の表面積及び流入エネルギーの経年変化を示すグラフである。
【図７】別の患者の動脈瘤の表面積及び流入エネルギーの経年変化を示すグラフである。
【図８】別の患者の動脈瘤の表面積及び流入エネルギーの経年変化を示すグラフである。
【図９】複数の患者のそれぞれについての「動脈瘤の表面積増加量と流入エネルギー（圧力エネルギーと運動エネルギーとの和）との関係」を示すグラフである。
【図１０】複数の患者のそれぞれについての「動脈瘤の表面積増加量と流入エネルギー（圧力エネルギー）との関係」を示すグラフである。
【図１１】本発明の実施形態に係る流体解析に使用される要素の斜視図である。
【図１２】本発明の実施形態に係る流体解析に使用される要素の斜視図である。
【図１３】本発明の実施形態に係る流体解析において採用されたＳＩＭＰＬＥ法による処理手順を示すフローチャートである。
【図１４】本発明の実施形態に係る流体解析に使用される「境界面に接する要素」の斜視図である。
【図１５】本発明の実施形態に係る流体解析に使用される「境界面に接する要素」の断面図である。
【図１６】本発明の実施形態に係る流体解析における処理手順を示すフローチャートである。
【発明を実施するための形態】
【００１５】
以下、本発明による「動脈瘤成長判定方法及び同方法に使用される装置」の実施形態について図面を参照しながら説明する。
【００１６】
＜画像取得＞
先ず、ＣＴスキャン（コンピュータ断層撮影）により「動脈瘤（本例では、脳動脈瘤）及びその周辺部」の３次元画像を取得する。図１は取得された３次元画像の一例を模式的に示す。この例において、動脈瘤１０は、複数の血管１１〜１５が集合している部分に発生している。血管１３は動脈瘤１０に流入する血液が流れる血管である。その他の血管（１１、１２、１４及び１５）は動脈瘤１０から流出する血液が流れる血管である。なお、図１に示した形状は一例であって、動脈瘤１０が発生する部位及び動脈瘤１０の周囲の血管の形状は様々であることは言うまでもない。また、上記画像は、ＣＴスキャンに限られることはなくＭＲＩ等を用いて取得されてもよい。
【００１７】
＜血液の流速及び圧力の取得＞
次に、動脈瘤１０の血液流入口１０ａにおける「血液の流速ｖin及び圧力Ｐin」を取得する。本実施形態において、これらの値は、後述する流体解析（非定常流体解析）により求められる。なお、係る解析法の基礎事項は良く知られている（例えば、「数値流体力学：Versteeg & Malalasekera著：森北出版」、「数値流体力学：越塚誠一著：培風館」及び「数値流体工学：荒川忠一著：東京大学出版会」等を参照。）。
【００１８】
この解析においては、動脈瘤及びその動脈瘤に接続した血管を「非圧縮性流体としての血液」が流れる構造体として扱い、その血管の内部を流れる血液に有限体積法を適用することにより前記血液の流速及び圧力を取得する。即ち、その非圧縮性流体の基礎方程式を数値的に解くために有限体積法が採用されている。有限体積法に用いられるメッシュとしては、血管（壁）付近には５面体要素が用いられ、血管中心部には４面体要素が用いられる。
【００１９】
加えて、図２に示したように、動脈瘤１０から流出する血液が流れる血管１１、１２、１４及び１５のそれぞれの端部１１ａ、１２ａ、１４ａ及び１５ａにおける「境界条件としての圧力」は「０」に設定される。なお、動脈瘤１０に流入する血液が流れる血管１３の端部１３ａの圧力は後述する解析により計算により求められる。
【００２０】
更に、動脈瘤１０に流入する血液が流れる血管１３の端部（流入口端）１３ａにおける「境界条件としての流速」は「所定流速」に設定される。この流速はベクトルであり、方向は端部１３ａと直交する方向（血管１３の中心軸に平行な方向）である。
【００２１】
この「所定流速（速度）」は次のようにして決定される。
（１）「太さ及び身体（例えば、頭部）における位置」が互いに異なる種々の血管について、その血管を流れる血液の流速を予め測定しデータ化しておく。
（２）上述したＣＴスキャンによって取得された「動脈瘤１０の画像」から、血管１３の太さ及び端部１３ａの身体における位置を特定する。
（３）特定された太さ及び位置に対応する（相関を有する）血管内を流れる血液の流速を「予め測定しておいた前記データ」に基づいて決定する。この決定された流速が前述した「所定流速」として採用される。
【００２２】
例えば、図２のグラフは、血管１３と同等の太さを有する血管であって、端部１３ａの近傍部位を通過する血管についての「一回の脈動（拍動）における血液の流速」の測定結果を示す。このような流速が境界条件として与えられる。
【００２３】
更に、この解析においては、血液表面（即ち、血管壁と接する部分）の流速は「０」に設定される。
【００２４】
このような境界条件を用いて流体解析（非定常流体解析）を行うことにより、動脈瘤１０の血液流入口（血液流出入断面部）１０ａにおける「血液の流速（ベクトル）ｖin及び圧力Ｐin等」が取得される。なお、血液流入口１０ａは、図１及び図２に示したように、動脈瘤１０に接続している血管と動脈瘤１０との境界平面（平面領域）を意味する。
【００２５】
＜動脈瘤へ流入する血液の流入エネルギーＥの取得＞
次に、図３に示したように、流速（ベクトル）ｖinのうち、動脈瘤１０に流入する血液の血液流入口（境界平面）１０ａに直交する成分（流速面直成分）ｖ’in（ベクトル）を求める。次に、血液流入口（境界平面）１０ａを通過して動脈瘤１０に流入する血液の体積流量Ｑinを下記（１）式に基づいて求める。（１）式においてＡinは「流速面直成分ｖ’inの血液」が流れ込む境界平面１０ａの面積である。
【数１】
2015073858.tif 000003
【００２６】
次いで、下記（２）式に基づいて、動脈瘤１０に流入する流入エネルギーＥを求める。（２）式において、ρは血液の密度である。（２）式の右辺第１項は血液流入口（境界平面）１０ａを通って動脈瘤１０に流入する血液の運動エネルギーであり、右辺第２項はその血液の圧力エネルギーである。
【数２】
2015073858.tif 000004
【００２７】
なお、実際の計算においては、有限体積法の各要素の血液流入口（境界平面）１０ａに沿う微小な面領域毎に流入エネルギーＥｉが計算され、その流入エネルギーＥｉが境界平面１０ａのうちの血液が動脈瘤１０に流入している領域全体に渡って積算されることにより、流入エネルギーＥが求められる。また、後述するように、流入エネルギーＥとして、動脈瘤１０に流入する血液の圧力エネルギー（Ｐin・Ｑin）を採用してもよい。
【００２８】
＜流入エネルギーＥと動脈瘤の成長との関係＞
次に、発明者は、流入エネルギーＥと動脈瘤（脳動脈瘤）１０の成長との関係について調査した。その結果を、図４乃至図８に示す。図４乃至図８は、各年における患者の「上記流入エネルギー及び動脈瘤の表面積からなるポイント」をプロットしたグラフである。但し、図４乃至図８における流入エネルギーは、運動エネルギーと圧力エネルギーとの和である。動脈瘤の表面積には血液流入口（境界平面）１０ａの面積は含まれない。
【００２９】
図４及び図５は、流入エネルギーＥが比較的大きい患者（Ｎ１，Ｎ２）についてのデータである。例えば、これらの患者の流入エネルギーＥは１．５（Ｎ・ｍ／ｓ）以上である。これらの患者については、動脈瘤の表面積が次第に増加していること（即ち、動脈瘤が成長していること）が判明した。
【００３０】
これに対し、図６乃至図８は、流入エネルギーＥが比較的小さい患者（Ｎ３、Ｎ４及びＮ５）についてのデータである。例えば、これらの患者の流入エネルギーＥは１．５（Ｎ・ｍ／ｓ）未満である。これらの患者については、動脈瘤の表面積は略一定であり、殆ど増加していないこと（即ち、動脈瘤が成長しておらず安定していること）が判明した。
【００３１】
図９は、動脈瘤の表面積増加量と流入エネルギーＥとの関係を示したグラフである。但し、図９における流入エネルギーＥは、前述した動脈瘤に流入する血液の「圧力エネルギーと運動エネルギーとの和」である。図１０は、動脈瘤の表面積増加量と流入エネルギーＥとの関係を示したグラフである。但し、図１０における流入エネルギーＥは、前述した動脈瘤に流入する血液の「圧力エネルギー」である。また、これらの図における表面積増加量は、１年（一定時間）当たりの「動脈瘤の表面積の増加量」である。
【００３２】
図９及び図１０から、流入エネルギーＥが大きいほど動脈瘤の表面積増加量が大きくなるという結論が得られた。従って、本実施形態においては、流入エネルギーＥに基づいて動脈瘤が成長するか否かを判定する。即ち、流入エネルギーＥが大きいほど、動脈瘤が成長する可能性が高いと判定する。これにより、動脈瘤が成長する可能性を精度よく推定することができる。
【００３３】
別の方法では、流入エネルギーＥが所定の閾値Ｅｔｈ（例えばＥｔｈ＝１．５（Ｎ・ｍ／ｓ））以上であるか否かを判定し、流入エネルギーＥが所定の閾値Ｅｔｈ以上であるとき、動脈瘤が成長すると判定する。一方、流入エネルギーＥが所定の閾値Ｅｔｈ未満であるとき、動脈瘤は安定していると判定する。これによっても、動脈瘤が成長するか否かを精度よく判定することができる。
【００３４】
＜血液の流速ｖin及び圧力Ｐinを取得するための流体解析＞
ところで、上述したように、上記（１）式及び（２）式にて流入エネルギーＥを算出するために用いられる「流速ｖin及び圧力Ｐin」は流体解析（非定常流体解析）に基づく数値計算をコンピュータにより行うことによって求められる。以下、この解析手法について述べる。但し、以下の手法は一例であり、本願発明がこれに限定されるものではない。
【００３５】
１．流体の理論式
この解析においては、血液を非圧縮性流体として扱う。更に、その非圧縮性流体は、基礎方程式としての「連続の式及び運動量保存の式」を満足するとの仮定に立つ。これらの式は周知である。なお、運動量保存の式においては、血液の粘性も考慮される。一方、温度分布は存在していないと考えてよいので、エネルギー保存の式は無視する。更に、血液の流れはレイノルズ数が低いので、層流として扱う。
【００３６】
上記連続の式は次の（１−１）式である。
【数３】
2015073858.tif 000005
【００３７】
上記運動量保存の式は、次の（１−２）、（１−３）及び（１−４）式である。運動量保存の式は「Navier-Stokesの式」とも称呼される。
【数４】
2015073858.tif 000006
【数５】
2015073858.tif 000007
【数６】
2015073858.tif 000008
【００３８】
ここで、
｛ｕ｝＝（ｕ，ｖ，ｗ）は速度（[ｍ／ｓ]）であり、
Ｐは圧力（[Ｐａ]＝[Ｎ／ｍ２]＝[Ｋｇ／（ｍ・ｓ２）]）であり、
ρは密度（[Ｋｇ／ｍ３]）であり、
μは粘性係数（[Ｋｇ／（ｍ・ｓ）]）である。
【００３９】
２．数値解法：有限体積法
２−１：有限体積法の理論式
上述したように、これらの基礎方程式を解くために、本例では有限体積法を採用した。なお、計算精度の確保及びメッシュの自動作成を目的として、血管近傍部位は５面体要素（プリズム要素）を用いるとともに血管中心部は４面体要素（テトラ要素）を用いてモデル化を行った。
【００４０】
より具体的に述べると、上記連続の式（１−１）を流体の領域内部の微小空間ＣＶで積分し、更に、Ｇａｕｓｓの発散定理を適用することにより、次の（２−１）式が得られる。なお、{ｎ}は図１１に示したように、微小空間の表面Ｓにおける法線ベクトルである。
【数７】
2015073858.tif 000009
【００４１】
同様に、上記運動量保存の式（１−２〜４）式からも、以下の（２−２〜４）式が得られる。
【数８】
2015073858.tif 000010
【数９】
2015073858.tif 000011
【数１０】
2015073858.tif 000012
【００４２】
２−２：離散化
本例では、数値計算のための離散化を行うため、速度及び圧力は各要素の中心位置で定義されるコロケート格子とする。
【００４３】
更に、本例では、要素（４面体要素及び５面体要素）毎に基礎方程式を体積積分することにより、要素毎の圧力と流速とを未知数とする非定常方程式（下記の（２−６）式）を得て、計算を行う。
【００４４】
具体的に述べると、要素の界面における流速（即ち、要素の界面における流量）Ｆｓを次の（２−５）式により定義する。ここで、ΔＡｓは界面の面積、{ｕｓ}は界面上の速度、{ｎｓ}は界面外向きの法線ベクトルである。
【数１１】
2015073858.tif 000013
【００４５】
すると、運動量保存の式（Navier-Stokesの式、（２−２〜４）式）は、次の（２−６）式に置換される。ここで、ΔＶは要素の体積、Ｐｓは界面上の圧力、[τｓ]は界面上の応力テンソルである。
【数１２】
2015073858.tif 000014
【００４６】
３．非定常流れに対する有限体積法
本例では血液を非圧縮性流体として扱っているので、連続の式である（１−１）式には時間微分項がないが、運動量保存の式（２−２〜４）式には時間微分項がある。時間微分項を扱う手法には、陽解法、Ｃｒａｎｋ−Ｎｉｃｏｌｓｏｎ法（中心差分法）、及び、多用される陰解法等があり、何れを採用してもよい。但し、本例では、コーディングが容易であり且つ連立一次方程式の解法が不要である陽解法を採用した。
【００４７】
具体的に述べると、置換された運動量保存の（２−６）式に陽解法を適用することにより、次の（３−１）式が得られる。この式において、右下の添え字「Ｐ」は着目要素を表し、右上の添え字「ｎ」は時刻を表す。右辺第２項の括弧内の１番目の項は対流項、２番目の項は圧力項、３番目の項は拡散項である。
【数１３】
2015073858.tif 000015
【００４８】
陽解法を適用する場合、時刻刻みΔｔは次の（３−２）式により示されるクーラン条件を満足する必要がある。これは、時刻刻みΔｔの間に許される流体の移動距離が１メッシュ未満であることに相当する。例えば、最大速度を８００[ｍｍ／ｓ]、メッシュサイズを０．２[ｍｍ]とすると、Δｔは０．２５[ｍｓ]以下でなければならない。
【数１４】
2015073858.tif 000016
【００４９】
４．拡散項の離散化
（３−１）式における拡散項における「界面上の応力テンソル[τｓ]と、図１２に示した界面を挟む要素間の中心間を結ぶベクトル｛γPA｝と、の内積は次の（４−１）式により計算することができる。なお、μは粘性係数である。
【数１５】
2015073858.tif 000017
【００５０】
ここで、ベクトル｛γPA｝の向きと界面の法線ベクトル｛ｎs｝の向きとが等しいと仮定すると、（３−１）式における拡散項は次の（４−２）式で表される。
【数１６】
2015073858.tif 000018
【００５１】
上記仮定は一般的には成立しないが、以下の理由により（４−２）式を採用する。
（１）拡散項は流体解析結果に与える影響が小さい。
（２）メッシュに分割した結果として生じる歪み及びアスペクト比が適正な範囲であれば、上記仮定は成立する。
（３）計算量を少なくすることができる。
【００５２】
（４−２）式を用いることにより、（３−１）式における拡散項は次の（４−３）式で表される。
【数１７】
2015073858.tif 000019
【００５３】
上記（４−３）式の係数は以下のとおりである。
【数１８】
2015073858.tif 000020
【数１９】
2015073858.tif 000021
【００５４】
５．対流項の離散化
（３−１）式における対流項は次式で与えられる。
【数２０】
2015073858.tif 000022
【００５５】
流束Ｆｓを計算するためには界面中心の速度（界面速度）ｕｆ＝｛ｕｓ｝・｛ｎｓ｝を計算する必要がある。本例では、ＲｈｉｅとＣｈｏｗの補間方法を適用した。これによれば、要素Ｐと要素Ａとの間の界面速度ｕｆは次式で与えられる。
【数２１】
2015073858.tif 000023
【００５６】
ここで、∇ｐ｜ｐは要素Ｐの中心における圧力の勾配であり、∇ｐ｜Ａは要素Ａの中心における圧力の勾配である。これらは、勾配再構築（再構成）法を用いることにより、接続する要素の圧力から計算することができる。
【００５７】
ａｐ及びａＡは、定常状態のＮａｖｉｅｒ−Ｓｔｏｋｅｓの式を離散化した際に、要素Ｐと要素Ａの速度に関する対角項となる値であり、（２−６）式と（４−２）式とから計算することができる。一例としてａｐの導出を以下に示す。
【００５８】
（２−６）式を定常状態の式に置き換えると次式が得られる。
【数２２】
2015073858.tif 000024
この式に（４−２）式を適用すると次式が得られる。
【数２３】
2015073858.tif 000025
【００５９】
更に、界面上の速度｛ｕｓ｝を、その界面に接続する要素の中心差分として与え、更に、両辺にρを乗じることにより次式が得られる。
【数２４】
2015073858.tif 000026
【００６０】
この式から、要素Ｐの速度｛ｕｐ｝に関する係数を抽出すると、次の（５−３）式が得られる。なお、（５−３）式は風上差分法ではなく中心差分法を適用している。
【数２５】
2015073858.tif 000027
以上により、｛ｕｆ｝が求められるので、（２−５）式から流束Ｆｓを計算することができる。
【００６１】
一方、（５−１）式の対流項の速度｛ｕｓ｝の計算には、精度向上のために中心差分法ではなく風上差分法を適用する。本例においては、風上差分法としてＴＶＤ（Ｔｏｔａｌ  Ｖａｒｉａｔｉｏｎ  Ｄｉｍｉｎｉｓｈｉｎｇ）法を用いる。このＴＶＤ法においては、要素Ｐと要素Ａとの界面におけるφｓを次式により与える。ここで、Ｕは風上側要素を表し、Ｄは風下側要素を表す。
【数２６】
2015073858.tif 000028
【００６２】
ここで、γｓは要素Ｐと要素Ａとの間のベクトル｛γPA｝を用いた次式で定義される。
【数２７】
2015073858.tif 000029
【００６３】
一方、（５−４）式の流速制限関数には以下に示す複数の関数の中の一つを採用すればよい。
【数２８】
2015073858.tif 000030
【００６４】
（５−５）式における∇φＰは、以下に述べるように、勾配再構築のための最小２乗法により求める。
【００６５】
要素Ｐが要素Ａ，Ｂ，Ｃ，Ｄに界面にて接しているとすると、次の式（（５−１３〜１６）式）が得られ、これらを行列表現すると（５−１７）式が得られる。
【数２９】
2015073858.tif 000031
【数３０】
2015073858.tif 000032
【００６６】
（５−１７）式の最小２乗解は次式により得られる。これにより得られる最小２乗解を∇φＰとして使用する。
【数３１】
2015073858.tif 000033
【００６７】
ところで、（５−４）式をそのまま用いると、界面速度が変数として現れる。従って、以下の手順により界面速度を要素中心速度に置き換える。
先ず、（５−１）式に（５−４）式を適用する。
【数３２】
2015073858.tif 000034
【００６８】
この結果、（５−１）式は次式のように表される。
【数３３】
2015073858.tif 000035
【００６９】
（５−２２）式の係数は以下のように導出することができる。先ず、（５−１９）式を速度の方向成分毎に記述すると次の３つの（５−２３、２４、２５）式が得られる。
【数３４】
2015073858.tif 000036
【００７０】
（５−２３，２４、２５）式を行列により表現すると、次の式が得られる。
【数３５】
2015073858.tif 000037
【００７１】
この式から、（５−２２）式の各係数は次の（５−２７）式及び（５−２８）式により表される。
【数３６】
2015073858.tif 000038
【数３７】
2015073858.tif 000039
【００７２】
６．圧力項の離散化
（３−１）式における圧力項は次の（６−１）式で与えられる。
【数３８】
2015073858.tif 000040
【００７３】
本例においては、この圧力項を界面上の圧力Ｐｓを使った次の（６−２）式により求める。この式は近似式である。
【数３９】
2015073858.tif 000041
【００７４】
ここで、界面上の圧力Ｐｓは、次の（６−３）式に示したように、その界面を共有する２つの要素の圧力の平均として近似的に求める。
【数４０】
2015073858.tif 000042
【００７５】
以上から、（６−１）式は次の（６−４〜６）式のように変形される。これらの式により圧力項が求められる。
【数４１】
2015073858.tif 000043
【００７６】
７．ＳＩＭＰＬＥ法による圧力の更新
上記（３−１）式を用いることにより時刻ｎ＋１の速度を計算することができる。しかしながら、圧力そのものは連続の式を用いて以下のように計算する。この圧力の計算には、ＳＩＭＰＬＥ（Semi-Implicit Method for Pressure-Linked Equations）を適用する。ＳＩＭＰＬＥ法は、以下に述べるように速度と圧力の補正値を用いる。
【００７７】
先ず、（３−１）式に陰解法を適用し、更に、（４−３）式、（５−２２）式及び（６−４）式を用いて書き換えると、次の（７−１）式が得られる。
【数４２】
2015073858.tif 000044
【００７８】
ここで、（３−１）式にて予測した速度を｛ｕＰ｝＊とすると、（７−１）式は次の（７−２）式に変形される。
【数４３】
2015073858.tif 000045
【００７９】
ここで、真の速度ｕ、予測した速度ｕ＊、速度の補正値ｕ’の間には下記の（７−３）式の関係があると仮定する。
同様に、真の圧力ｐ、予測した圧力ｐ＊、圧力の補正値ｐ’の間には下記の（７−４）式の関係があると仮定する。
【数４４】
2015073858.tif 000046
【００８０】
（７−１）式から（７−２）式を差し引き、（７−３）式及び（７−４）式を適用すると次の（７−５）式が得られる。
【数４５】
2015073858.tif 000047
【００８１】
（７−５）式中の隣接要素速度に関する和を無視すると次の式が得られる。
【数４６】
2015073858.tif 000048
【００８２】
（７−８）式の逆行列は、対角行列（対角行列である元の行列）の逆行列である。従って、（７−８）式の逆行列は、元の行列の対角項を逆数にすることにより求めることができる（即ち、ベクトルの成分毎に逆数を乗じればよい。）。
なお、以降において出現する記号の意味を次の（７−９）式により表しておく。
【数４７】
2015073858.tif 000049
【００８３】
上記（７−７）式は、「速度の補正値は、隣接する要素の圧力の補正値と一致する。」ことを示している。
【００８４】
更に、（７−７）式の関係を連続の式に適用することにより、圧力の補正値に関する式を導出する。
先ず、要素Ｐに関する連続の式である（２−１）式は次の式により表される。
【数４８】
2015073858.tif 000050
【００８５】
ここで、次式に示したように、界面中心速度｛ｕｓ｝ｎ＋１は隣接する要素の速度の平均で表されるものと仮定する。
【数４９】
2015073858.tif 000051
【数５０】
2015073858.tif 000052
【００８６】
ここで、閉空間の表面積と法線ベクトルとの間には次の関係式が成立する。
【数５１】
2015073858.tif 000053
【００８７】
これを（７−１２）式に適用すると次の（７−１４）式が得られる。
【数５２】
2015073858.tif 000054
【００８８】
これに（７−７）式を適用すると、次の（７−１５）式が得られる。
【数５３】
2015073858.tif 000055
【００８９】
（７−１５）式を整理すると次の（７−１６）式が得られる。（７−１６）式の係数は（７−１７）式及び（７−１８）式により表される。
【数５４】
2015073858.tif 000056
【数５５】
2015073858.tif 000057
【数５６】
2015073858.tif 000058
【００９０】
（７−１６）式を各要素について連立させれば、圧力の補正値ｐ’に関する連立一次方程式が得られる。これを解くことにより圧力の補正値を得ることができる。この連立一次方程式はポアソン方程式の解法に相当していると考えられる。
更に、（７−１６）式を解くことにより得られた圧力の補正値ｐ’を（７−７）式に代入することにより、要素中心の速度の補正値｛ｕＰ｝’を得ることができる。
【００９１】
以上により得られた補正値を用いて圧力と速度を更新する。但し、圧力の補正値の発散を防ぐためにここでは緩和法を適用する。緩和係数αｐを（０＜αｐ＜１）の値とすると、圧力の更新式は次の（７−１９）式のように表現される。
【数５７】
2015073858.tif 000059
【００９２】
速度も緩和係数を用いて更新を行う。但し、次の（７−２０）式を用いる。
ここで、｛ｕＰ｝oldは、ＳＩＭＰＬＥ法における一回前の反復計算により求めた値である。
【数５８】
2015073858.tif 000060
【００９３】
上記（７−２０）式により更新した速度を下記の（７−２１）式に適用して流束Ｆｓnewを計算する。そして、（７−１０）式の連続の式が満足されなければ、（７−１６）式を用いた圧力補正値の計算を繰り返し行う。
【数５９】
2015073858.tif 000061
【００９４】
以上に説明したＳＩＭＰＬＥ法の計算手順を参考として図１３のフローチャートにより示した。流速が連続の式を満たしているか否かの判定は、数値誤差を考慮して最大の流速の相対値を基準に判定している（Ｓ４を参照。）。なお、｛ｕＬ｝は流入速度である。
【００９５】
８．境界条件
本例においては、境界条件として、次の３種類の条件を設定できる。
（１）速度の時間変化
（２）圧力の時間変化
（３）壁面
なお、壁面は可動であるが、速度に比べ壁面の移動速度は微小であるので、無視する。よって、構造側との連成は圧力の釣り合いを満足することにより実現する。また、固体の壁での速度成分に対しては、すべりなしの条件（non-slip condition) ｛ｕ｝＝０が適切
な条件であり、これを本例にも適用した。
【００９６】
（１）拡散項の境界条件
（１Ａ）速度境界条件
前述した「４．拡散項の離散化」において、要素Ｐを速度境界面に接続する要素とする。そして、要素Ｐに接続する要素Ａの代わりに、図１４に示したように、要素Ｐの中心から境界面に下した垂線の足の位置をＬとすると、（４−２）式は次の（８−１）式となる。｛ｕＬ｝は境界面の速度である。この（８−１）式を（３−２）式に適用することにより、速度境界条件を取り込むことができる。
【数６０】
2015073858.tif 000062
【００９７】
（１Ｂ）圧力境界条件
圧力境界条件が設定された場合、境界面の速度｛ｕＬ｝は、要素の中心の速度｛ｕｐ｝n+1と等しいと考える。よって圧力境界条件として次の（８−２）式が得られる。
【数６１】
2015073858.tif 000063
【００９８】
（１Ｃ）壁面境界条件
上述のすべりなしの条件（non-slip condition)を適用すると、｛ｕＬ｝＝０である。よって、（８−１）式は次の（８−３）式となる。
【数６２】
2015073858.tif 000064
【００９９】
（２）対流項の境界条件
（２Ａ）速度境界条件
速度境界の速度を｛ｕＬ｝、速度境界の要素から見た境界面の法線ベクトルを｛ｎＬ｝とするとき、速度境界の界面速度ｕｆは（５−２）式に代わり、次の（８−４）式により表される。
【数６３】
2015073858.tif 000065
【０１００】
速度境界における流束Ｆｓは、この速度ｕｆを用いて計算する。更に、（３−２）式における界面の速度｛ｕＳ｝n+1は、速度境界の速度｛ｕＬ｝をそのまま利用することができる。
【０１０１】
ところで、境界要素Ｐのａｐは境界要素Ｐに接続する要素との界面速度計算及びＳＩＭＰＬＥ法での連立一次方程式の作成に必要である。境界要素Ｐのａｐも上記（５−３）式を用いて求める。但し、境界面の速度｛ｕＳ｝が既知（｛ｕＳ｝＝｛ｕＬ｝）であるから、境界要素Ｐのａｐを求める場合の（５−３）式は次の（８−５）式へと変更される（（８−５）式の変数については図１５を参照）。
【数６４】
2015073858.tif 000066
【０１０２】
（２Ｂ）圧力境界条件
圧力境界条件が設定された場合、境界面の速度｛ｕＬ｝は要素中心の速度｛ｕＰ｝と等しいと考えるので、流束ＦＬを計算する際のｕｆは次の（８−６）式から計算する。
【数６５】
2015073858.tif 000067
【０１０３】
（３−２）式における界面の速度｛ｕＳ｝n+1には要素中心の速度｛ｕＰ｝n+1をそのまま利用すればよい。
更に、境界面での速度｛ｕＳ｝＝｛ｕＰ｝であるから、（５−３）式は次の（８−７）式に示すようになる。ＦＬは境界面の流束である。
【数６６】
2015073858.tif 000068
【０１０４】
（２Ｃ）壁面境界条件
上述のすべりなしの条件（non-slip condition)を適用すると、｛ｕＬ｝＝０である。よって、（８−１）式は次の（８−８）式となる。
【数６７】
2015073858.tif 000069
【０１０５】
（３）圧力項の境界条件
境界条件が設定されている場合（即ち、隣接要素が３個以上ない場合）、（６−１）式で必要な圧力勾配∇ｐの再構築のために最小２乗法を適用できない。よって、この場合、圧力勾配はゼロ（∇ｐ＝０）とする。
【０１０６】
（４）ＳＩＭＰＬＥ法の取り扱い
（４Ａ）速度境界条件
速度境界の流束（既知）をＦinとすると、（７−４）式〜（７−７）式は次のとおりとなる。ここで、inner
surfacesとは速度境界以外の界面を意味する。
【数６８】
2015073858.tif 000070
【０１０７】
また、速度境界では圧力は変化しないと仮定すればよい。よって、圧力勾配を（６−４）式にて計算する際の界面圧力ＰＳは、次の（８−１４）式により計算される。
【数６９】
2015073858.tif 000071
【０１０８】
（４Ｂ）圧力境界条件
圧力境界条件が設定された場合、境界上の要素中心の圧力ＰＰは境界上の圧力ＰＬと等しいと考える。
（４Ｃ）壁面境界条件
壁面への流入出はないので、（７−７）式の和から壁面を除外する。また、壁面境界では圧力は変化しないと仮定すればよい。よって、圧力勾配を（６−４）式により計算する際の界面圧力ＰＳはＰＳ＝ＰＰとして計算する。
【０１０９】
＜実装＞
一般に、このような有限体積法を用いた解析（シミュレーション）は、入力部、記憶部、演算部及び出力部をハードウェアとして有する汎用又は専用コンピュータにより実施することができる。即ち、入力部を用いて「流体（血液）の流れる構造物」を、複数の要素に分割することによりモデルを構築し、そのモデルが記憶部に記憶され、演算部がそのモデルを用いて流体解析を実施する。出力部はその解析結果を出力する。なお、解析の前段の処理としてプリプロセッサを用いてもよく、解析の後段の処理としてポストプロセッサを用いてもよい。
【０１１０】
＜処理手順＞
図１６は以上に述べた解析に従って計算を行う処理の流れを示す概略フローチャートである。以下、各処理について簡単に説明する。
ステップ１０５：演算部は、血液の速度（流速）｛ｕ｝及び圧力Ｐの初期値を設定する。例えば、速度｛ｕ｝の初期値は、上述したように決定された「所定流速」の一回の脈動（拍動）における最大値に設定されてもよく、「所定流速」の一回の脈動における平均値に設定されてもよい。圧力Ｐの初期値は、所定の一定値（例えば、「０」又は患者の血圧に応じた「０以外の値」）に設定されてもよく、血管の部位毎の圧力を予め推定又は取得できる場合には、それらの部位毎の圧力に設定されてもよい。
【０１１１】
ステップ１１０：演算部は、時刻ｔをΔｔだけ進める。
ステップ１１５：演算部は、上記（４−１）式〜（４−５）式に従って拡散項の係数を演算する。
ステップ１２０：演算部は、上記（５−１）式〜（５−２８）式に従って対流項の係数を演算する。
ステップ１２５：演算部は、上記（６−１）式〜（６−６）式に従って圧力項の係数を演算する。
ステップ１３０：演算部は、時刻ｎ＋１（即ち、ｔ＋Δｔ）における血液の速度（流速）を上記（７−１）式〜（７−９）式に従って演算する。
【０１１２】
演算部は、ステップ１３５乃至ステップ１５０の処理により、図１３を参照しながら詳述したＳＩＭＰＬＥ法に従って所定の演算を行う。ここでは、図１３の各処理をより簡素化して記載した。
ステップ１３５：演算部は、上記（７−１０）式〜（７−１８）式に従って連続の式の係数行列及び右辺ベクトルを作成する（特に、（７−１７）式及び（７−１８）式を参照。）。
ステップ１４０：演算部は、上記（７−１６）式の連立一次方程式を解いて、圧力補正値を求める。
ステップ１４５：演算部は、上記（７−１９）式〜（７−２１）式に従って速度及び圧力を補正する。
ステップ１５０：演算部は、収束しているか否かを判定する（図１３のＳ２〜Ｓ４を参照。）。演算部は、収束していなければステップ１３５に戻る。演算部は、収束している場合、ステップ１５５に進む。
【０１１３】
ステップ１５５：演算部は、演算の終了時刻が到来したか否かを判定する。終了時刻は任意の時刻に設定される。演算部は、終了時刻が到来していなければ、ステップ１１０に戻る。演算部は、終了時刻が到来していれば、ステップ１９５に進んで演算を終了する。これらの演算の結果、血液流入口１０ａにおける「血液の流速ｖin及び圧力Ｐin」が演算され、流入エネルギーＥが演算される。
【０１１４】
以上、説明したように、本発明の実施形態によれば、動脈瘤に流入するエネルギーに基づいて動脈瘤が成長するか否かを精度よく判定することができる。。
【０１１５】
本発明は上記実施形態に限定されることはなく、本発明の範囲内において種々の変形例を採用することができる。例えば、流入エネルギーを求めるために必要な「血液の流速（ベクトル）ｖin及び圧力Ｐin」は上記実施形態とは別の解析手法を用いて計算されてもよい。
【符号の説明】
【０１１６】
１０…動脈瘤
１０ａ…血液流入口（境界平面）
１１−１５…血管
【図１】
2015073858.tif 000072
【図２】
2015073858.tif 000073
【図３】
2015073858.tif 000074
【図４】
2015073858.tif 000075
【図５】
2015073858.tif 000076
【図６】
2015073858.tif 000077
【図７】
2015073858.tif 000078
【図８】
2015073858.tif 000079
【図９】
2015073858.tif 000080
【図１０】
2015073858.tif 000081
【図１１】
2015073858.tif 000082
【図１２】
2015073858.tif 000083
【図１３】
2015073858.tif 000084
【図１４】
2015073858.tif 000085
【図１５】
2015073858.tif 000086
【図１６】
2015073858.tif 000087