脳動脈瘤の血流モデリング

メッシュ参考文献
定常流

脳動脈瘤は血管壁の内層の衰えによって引き起こされる異常な拡張である．動脈瘤は脳内で起こる．最大のリスクは，死亡率がほぼ40%に及ぶ破裂である．正確な診断は年齢や動脈瘤の位置等さまざまな要因によって異なる．

この応用例の目的は，CTスキャンで得られた患者特定の動脈瘤における血流を調べ，さまざまな血流力学パラメータに基づいて最も破裂が起こりそうな場所を予測するというものである．このモデルは血液の非ニュートン流体流れ特性を考慮に入れる．拍動流データは中大脳動脈におけるドプラ超音波法で試験的に得られる．モデルと形状の複雑性により，このモデルは計算時間とメモリの両方の点で高価なものになる．

この応用例は2つの部分に分かれている．最初の部分では一人の患者特有の動脈瘤における定常流を調べる．次に動脈瘤の拍動流に焦点を当て，破裂場所の予測を行う．この応用例は[Brunátová, 2021]に基づいている．

単位

Wolfram言語における単位の標準系はSI単位である．血流シミュレーションでは，センチ，グラム，秒を使うCGS単位系が主に使われ，質量密度はグラム毎立方センチメートルの単位になる．圧力の基本単位はダイン毎平方センチメートル[]であり1 =0.1` を使って標準SI単位に変換することができる．動的粘度の基本単位はポアズで1 =0.1` である．

Wolfram言語における有限要素法のコンテキストでの単位についての詳細はこちらに記載されている．

有限要素法パッケージをロードし$HistoryLengthを設定する：

メッシュ

まずメッシュを生成する．この場合メッシュはすでに存在し，ここから利用することができるので最初から生成する必要はない．しかし，元のメッシュは非常に大きい．計算時間の観点から，計算を加速する縮小されたメッシュが提供されている．しかし後に分かるように，これが結果の正確さにマイナスの影響を及ぼすことがある．より正確な解を得るためには，もとのメッシュを使うことができる．

境界メッシュをインポートする：
境界メッシュを可視化する：

境界のさまざまな点における境界条件を規定するために，境界の部分を区別しなければならない．これは要素マーカーを使って行える．まず境界要素から始める．

境界のさまざまな部分における境界条件を規定するために，境界の部分が区別できなければならない．これは要素マーカーを使って行える．要素のタイプに従って，メッシュ要素マーカー，境界要素マーカー，点要素マーカーの3種類の要素マーカーがある．

すべての境界メッシュマーカーをリストする：
境界メッシュを異なる要素マーカーに対応する色で表示する：

左側の絵に示された境界の赤い部分は流入口を表す．これは血液が動脈瘤に流れ込む部分である．右の絵の境界の緑と青の部分は流出口を示す．これは血液が動脈瘤から流出する部分である．境界の白い部分は壁を表す．

ノイマン値境界マーカーを設定する：

境界要素マーカーは後でノイマン境界条件を規定する際に使う．しかしディリクレ境界条件は点要素で規定されるため，点要素についてこの手法を繰り返すことが必要になる．

すべての点要素マーカーをリストする：
境界の点要素を異なる点要素マーカーに対応する色で表示する：
ディリクレ条件境界マーカーを設定する：

最後に圧力に対する基準点を選ぶ．圧力について条件が与えられていない場合，解は一意にはならない．この問題を解くために，圧力の基準値を規定する点を選ぶ．

圧力の基準点を選び，その位置をメッシュに表示する：

ディリクレ境界条件は点要素で規定されるが圧力基準点は点要素ではない可能性がある．ゆえに最近傍の点要素を求める必要がある．

要素マーカーが1に等しい境界点要素をすべてリストする：
pointElementsの座標を得る：
最近傍点要素を求める：
新しい圧力基準点の位置を示す：

最後に完全な要素メッシュを生成する．

境界メッシュから完全メッシュを生成する：

定常流

まず簡単に済ませるために，定常流の場合を解く．直接時間依存の拍動流の場合を解くよりもこの方が簡単である．2つ目のステップとしてWolfram言語の今後のバージョンでは，拍動流の場合が加わる．これで利点が増える．この概念は簡単な設定で説明できる．

定常流シミュレーションの目的は，動脈瘤壁のせん断応力を計算することである．困難なのは，流体流れモデルにおける高非線形性である．高非線形のPDEが解けないのは珍しいことではない．一般に2種類の方法がある．一つは例えばParametricNDSolveValueからParametricFunctionを作成し，その非線形性を反復的に増加させるというものがある．このアプローチは超弾性材料モデルを収束させるために一般的に使われる．2つ目のアプローチはここで使うものであるが，定常モデルを時間依存にして時間パラメータ t を使って非線形性を上昇させ，定常状態に達するまで統合するというものである．したがって方程式は定常流の計算に使うものであっても，時間依存として設定する．このことで休止状態から始めて，時間の関数 t としてゆっくりと流入速度を上げていくことができる．NDSolve系の関数は適応ステップ統合法が実装されており，NDSolveは自動的に時間刻み幅を増加・減少させることができるので非常に使いやすい．

境界条件

境界には，壁面，流出口，流入口の3つの部分がある．これは指定しなければならない3つの境界条件に対応している．さらに基準圧力の条件も規定する必要がある．

壁面境界条件

ここではいわゆる滑りなし条件は壁面に規定する．つまり，壁面の速度は0でなければならない．流体は境界にくっついたままで滑らない，という理由でこの条件は滑りなしと呼ばれるのである．

境界で速度に対する滑りなし条件を定義する：

流出条件

流出口では血流を妨げるものがあってはならない．それゆえ流出口で牽引力0を指定する．この境界条件は自由流出境界条件とも呼ばれる．ノイマン0境界条件がデフォルトの境界条件なので，境界で何も指定されていないときは何も指定できないし，それで構わない．

流出口で自由流出条件を指定する：

基準圧力

次に圧力の境界条件を規定する．この境界条件は解を一意にするために必要である．

圧力基準点における圧力の条件を定義する：

流入条件

ここでは流入条件を定義するが，これはいくつかの部分からなる．

流入プロファイル

最初の部分は単位放物プロファイルである．これは流入の一般的なプロファイルを定義するのに使われ放物形である．放物形は中心に一つの値があり，隣接壁の滑りなし条件にマッチするように境界ではゼロに等しい．

これは管内のポアズイユ流れの解析解として与えられる．これは両端で異なる圧力を持つ円筒管の層流である．圧力差により管内の流体を流す圧力勾配が発生する．

流入プロファイルの式は以下で与えられる：

ここで r は管（流入口）の半径，は流入境界の中心点である．流入速度に対するこの式は，壁面に滑りなし条件があるときだけ有効である．異なる境界条件が使われると，この式は異なる形式になる．

この式を使うためには半径と中心を知る必要がある．ここではその値はメッシュで与えられる．

与えられた値から流入口の中心と半径を定義する：

ここで流入プロファイルが計算できる．プロファイルの定義に現れるベクトル(x-sx,y-sy,z-sz)は流入口にないかもしれない．結果のプロファイルを混乱させる可能性のある小さい法線成分がある可能性がある．この法線成分の可能性をベクトルから除去する．こうすることで流入口に接した，プロファイルの計算に使えるベクトルが得られる．このためには単位境界ベクトルが必要である．ここでも与えられたベクトルを使うがBoundaryUnitNormal関数という便利なツールが使える．

ベクトルから法線成分を除去することでベクトルの接線成分を計算する：
流入プロファイルを設定する：
流入減衰

流入減衰関数は，上で説明した通り最初に流入を減衰させるために使う．減衰関数の実際の数式は関数が0から始まる限り任意であり，1に達するまで増加しその後は一定になる．

減衰関数を定義する：
区間[0, 4]上の減衰関数をプロットする：
流入スケーリング

周期平均の断面平均速度を持つ流入をスケールする．名前の通りこれは流入口の断面および1心周期の平均速度である．その値は48 である．

周期平均の断面平均速度を定義する：
完全流入プロファイル

これで完全な流入関数が書ける．

流入速度関数を定義する：

この関数を使って，流入口の境界条件がついに定義できる．

流入口でディリクレ条件を定義する：
すべてのディリクレ境界条件を1つの変数にラップする：

初期条件

すべての境界条件が規定されたが，時間依存問題を解きたいので初期条件も与えなければならない．計算の最初では流体が休止した状態にしたいので，初期条件として速度も圧力もゼロを選ぶ．

初期条件を定義する：

モデルの設定と解

まず変数を定義する．

変数を定義する：

次に材料パラメータを定義する．人間の血液の特性を記述するには，非ニュートンCarreauモデルが適している．これを使うと，せん断速度が上昇したときの粘度の減少とニュートン領域を捉えることができる．

材料パラメータを定義する：

方程式はFluidFlowPDEComponentを使って生成することができる．

方程式を生成する：

これで方程式が解ける．これには有限要素法を使って方程式を空間領域に離散化し，線の方法を使って時間に対処する．また流体力学でよく使われるTaylor-Hood要素を使う．AccuracyGoalとPrecisionGoalを減少させ，次数の少ないIDA法を使って計算を高速にするが，正確性が減る．

さらに，計算が終了するときのためだけの解を保存するので十分である．このためNDSolveValueでを使う．これにより返されたInterpolatingFunctionオブジェクトを保存するのに必要なメモリが格段に少なくなる．

計算が終る時間を定義する：

シミュレーションは終了するまでに時間がかかる．

方程式を解き，進行状況を監視し，かかる時間を測定する：

メッシュは多数の要素を持った3次元であり数学モデルは非常に非線形なので，長い計算時間が予測される．これらの要因で評価時間が長くなる．

結果

これで方程式が解けたので，解を見て追加のパラメータを計算する．

速度

速度プロットを使ってどのように血液が流れるかを可視化する．

VectorPlot3Dを使って定常解をプロットする：

血液が流入口から動脈瘤に流れ，流出口から出るまで動脈瘤の中を回っている．この挙動はストリームプロットで見ることができる．ベクトルプロットでは，血液が流入口から動脈瘤に流れ込み流出口から流れ出るのが分かる．これはストリームプロットでは見られない想定通りの挙動である．

StreamPlot3Dを使って流線をプロットする：

もっと詳細を見るために，2つのスライスを作成し速度のノルムをプロットする．このスライスは法線ベクトルと点で定義される．これは動脈瘤をスライスするユニークな平面を定義する．

最初のスライスについて点と法線ベクトルを定義する：
2つ目のスライスについて点と法線ベクトルを定義する：
両方のスライスをプロットする：

両方のスライスを定義したので，これらのスライスに対して速度のノルムがどのように見えるかをプロットすることができる．

最初のスライスの速度のノルムをプロットする：
2つ目のスライスの速度のノルムをプロットする：

スライスの大きい血管部の中心の速度は非常に小さく速度のピークが境界付近で観察できる．

粘度

ここで速度のノルムをプロットしたのと同じ方法で見かけ粘度をプロットする．これにより，血液の非ニュートン特性が調べられる．しかし，粘度に直接アクセスできないので，まずこれを計算する必要がある．ひずみ速度テンソルから始めるが，そのためには速度勾配を計算する必要がある．

速度勾配を計算する：
ひずみ速度テンソルを計算する：

これでひずみ速度テンソルが得られたので，せん断速度を計算して見かけ粘度を計算することができる．

せん断速度を計算する：
見かけ粘度を計算する：
最初のスライスに見かけ粘度をプロットする：
2つ目のスライスに見かけ粘度をプロットする：

プロットは粘度が0.04` から0.1` まで変化するときを示しているので，血液の非ニュートン性を表す．せん断速度がゼロのときの粘度が0.56` であることを考慮に入れると，粘度が非常に減少している．ニュートン流体では両方のプロットが一定の粘度を示している．

壁面せん断応力

最後に定常解について行うのは，定常解から壁面せん断応力を計算することである．壁面せん断応力は，重要な血流力学パラメータである．これは壁面に作用する牽引力の接線部分を表し，破裂のサインの一つと考えられている．

壁面せん断応力の式は以下である：

ここで は応力テンソル， は外側の単位境界法線ベクトルである．最初の部分 は壁面に作用する牽引力を表し，2つ目の部分 は牽引の法線成分である．牽引から法線成分を引くと接線成分が残り，これを壁面せん断応力と言う．

壁面せん断応力について2つのことに気を付けなければならない．まず，ソースの中にはベクトル の法線として壁面せん断応力を定義するものもあれば，ベクトル自体として定義するものもある．

2つ目は，壁面せん断応力はメッシュ微調整に敏感ということである．この応用例では計算流体力学で発生する複雑な問題を解くのに関連した概念やワークフローを例示したいため，正確性よりも性能に焦点を当てている．このために縮小したメッシュを使い，壁面せん断応力について得られた結果は実際ほど正確ではないのである．問題の実践で使えるようなより正確な結果のためには，もとのメッシュ，理想的にはより細かいメッシュで問題を解く必要がある．

応力テンソルと単位境界法線ベクトルを設定して壁面せん断応力を計算しなければならない．応力テンソルは以下で与えられる：

ここで は圧力， は単位行列， は見かけ粘度， はひずみ速度テンソルである．前のセクションでこのほとんどは計算済みなので，応力テンソルと単位境界法線を設定する．

応力テンソル σ を定義する：
境界単位法線ベクトルを得る：

これで壁面せん断応力を計算しプロットすることができる．

壁に作用する牽引力を計算する：
壁面せん断応力を計算する：
壁面せん断応力をプロットする：

壁面せん断応力は，動脈瘤の破裂がどこで起こりやすいかを予測することにおいて重要な役割を果たす．破裂場所の予測は拍動流の場合に行うことができる．これはWolfram言語の将来のバージョンで加わる．

参考文献

1.  Brunátová (formerly Trdlicová), Jana. (2021). Blood Flow Modeling in Cerebral Aneurysm. Master's thesis, Charles University, Faculty of Mathematics and Physics, online available here.