平成 6 年度 卒業論文 狭窄部を有する血管内の血流の有限要素解析 高知工科大学工学部知能機械システム工学科知能流体力学研究室 清水昌彦
目次 第 章緒言 - 本研究を行う背景と目的 - 血液の性質 -3 数値計算 - 有限要素法の概要 第 章基礎方程式 - 支配方程式 -- 連続の式 5 -- コーシーの運動方程式 6 --3 血液の構成方程式 6 - 無次元化 7 第 3 章解析手法 3- 有限要素解析 8 3- 手順 8 3-- 時間方向の離散化 8 3-- 流速修正法 8 3--3 重み付き残差法 9 3-- 連立方程式の解法 3-3 プログラムの流れ 第 章計算条件と要素分割 - 計算条件 - 要素分割 第 5 章結果と考察 第 6 章結言 参考文献 謝辞 3
第 章緒言 - 本研究を行う背景と目的現在, 死亡原因の中の多くには心疾患や脳血管疾患がある. これらの病気の発生メカニズムは血液流動と血管に深く関わっており, 疾患と血液流動との関係を解明することが望まれている. しかし, 実際のヒトの血管を用いて実験することは困難であり, 流体力学の豊富な知識が必要とされるため, 医学分野のみで血液流動を解析するのは不可能である. よって, 流体力学の立場から考えることが必要である. そこで, 近年バイオメカにクスの分野では人体実験をする必要のないコンピュータシミュレーションによる解析は今まで行われてきたが, 血液をニュートン流体として解析されたものがほとんどで, 血液の非ニュートン性を考慮したものは少ない. そこで本研究では血液の非ニュートン性を考慮して血液流動の解析を行う. 解析対象として半球のついた円管内流れとし, 血管狭窄の状態に近づけた. この状況下でのニュートン流体, 非ニュートン流体の流動を解析し比較する. - 血液の性質血液は酸素や二酸化炭素の運搬や, 栄養物や老廃物の運搬などヒトにとって重要な役割を果たす. ヒトの体内を循環する血液は約 5% の細胞成分と約 55% の液体成分から構成される. 細胞成分には赤血球, 白血球, 血小板があり, 液体成分は血漿である. 血液に含まれる赤血球の容積比率は約 5% と非常に大きい. よって血液流動には赤血球の変形が大きく関わっている. 血液は水のようなニュートン流体ではなく, 血液に与えるせん断速度を大きくすれば見かけ粘度が小さくなる非ニュートン流体である. さらに血液の通り道である血管は粘弾性管である. 血液の非ニュートン性を表す式にキャッソンの式がある. キャッソンは印刷インクの流動挙動についてずり速度とずり応力の関係を円錐 - 平板型粘度計で測定して式を導いたが, この式が血液流動にも使えることがわかった. -3 数値計算数値計算とは解析的に解くことが不可能な自然現象などをモデル化し数値的に近似解を求め, 自然現象などを可視化することができる方法である. 対象によって数値計算の方法も使い分けなければならない. 本研究においては血液流動の解析を行うため, 実際の血管形状や, 疾患部を考慮し有限要素法を使うことにする. 流動の解析を行う際, 境界で変化が激しいことから有限要素法は計算領域に非構造格子を用いるため境界付近で格子を小さくし, 精度よく解析することが可能である.
- 有限要素法の概要 ある微分方程式 f がある. 近似解を とすると f ( cons. > ) n (.) r f (.) となり, r は残差である. これにある関数 をかけて積分したものを にしなければならない. よって, r (.3) となる. 関数 φ は重み関数といい, 重み関数によって得られる方程式を重み付き残差方程式という. 有限要素法では各要素を独立したものとして扱うので各要素によって近似関数を定義しなければならない. 各要素の近似関数は接点値を元に直線で近似する. この近似関数を要素補間関数という. 図. のように長さ の要素に着目し, 接点値を満たすように要素補間関数を定義すると b (.) ここで, は形状関数という. b b 図. 要素補間関数
3 また, 重み関数についても同様に近似関数を導入すると b (.5) となり, 式 (.),(.5) を式 (.3) に代入すると b b b (.6) 式 (.6) を積分すると b b b (.7) ここで重み関数を, b とすると b N (.8) また, b とすると b b N (.9) ここで N, b N はそれぞれ, b によって求められたことを表すために付け加えた. 式 (.8), (.9) をまとめると b b N N (.)
この式 (.) を有限要素方程式といい, この方程式は着目した要素内において接点の未知関数に関する量, b で離散化した方程式系となっている. この式を領域内すべての要素において求め, 有限要素方程式を領域全体で重ね合わせることで領域の支配方程式を単純な連立方程式に変換し解く. このように微分方程式を要素ごとに近似方程式にし, 領域全体で重ね合わせ, 支配方程式を近似した連立方程式を解くというのが有限要素法の原理である.
5 第 章基礎方程式 ここではまず主な記号について説明する. e : 変形速度テンソル (, ) 成分 e : 代表長さ : 代表速度 P : 圧力 : 降伏応力 S η : 時間 : キャッソン粘度 γ& : ずり速度 ( e ) ρ : 密度 τ : 偏差応力テンソル Π : e - 支配方程式 本研究において用いられている支配方程式は, 連続の式, コーシーの運動方程式およびキャッソン流体の構成方程式であり, 式の表記法として総和規約を用いる. また, 添字, は整数で ~3 の範囲を持つ.,, 3 のとき,,z 方向の方程式を表す. -- 連続の式連続の式は質量保存の式で流体が非圧縮のとき次式で表される. ここで, は座標成分であり,,, 3 z である. は速度の 成分である. (.)
6 -- コーシーの運動方程式コーシーの運動方程式は運動量保存の式であり, 領域に外部からなされる力がない場合, 次式のようになる. P τ ρ (.) --3 血液の構成方程式 ( ) τ S & γ (.3) η 上式 (.3) は血液の構成方程式としてよく用いられるキャッソンの式である. この式を τ について整理すると次式のようになる. S S η τ η Π Π e (.) しかし, Π のとき, τ が微分不可能になるという特性があるので修正キャッソンの式を導入し Π のときでも微分可能にするためにδ を加えた. S S η τ η e Π δ Π δ (.5) δ はずり速度が低い領域での血液粘度にあうように定めた. また, 血液の構成方程式の降 伏応力がないとき, つまり S のときはニュートン流体の構成方程式である.
7 - 無次元化 - の支配方程式を無次元化する. (.6) (.7) S S η (.8) P P ρ (.9) (.) e (.) Π Π (.) ここで代表長さ には血管の直径, 代表速度 には流入の平均流速を用いた. は無次元量を表している. 上式 (.6) から (.) を式 (.5) に代入すると e S S S S τ Π Π Π Π τ η η δ δ η δ δ (.3) 次にコーシーの運動方程式 (.) に式 (.6),(.7),(.9),(.) および (.3) を代入すると P τ η ρ ρ 上式を両辺 ρ で割る, すなわち, ρ をかけると Re P P τ τ ρ η (.) となる.(Re: レイノルズ数 )
8 第 3 章解析手法 3- 有限要素解析本研究では重み付き残差法によって数値解析を行う. 重み付き残差法とは微分方程式の残差と重み関数の内積を零に近づけていく方法である. ここではその実際の手順を述べる. 3- 手順 3-- 時間方向の離散化コーシーの運動方程式に対して時間方向の離散化 new new P Δ τ Re (3.) new (3.) 3-- 流速修正法流速修正法とは流速の予測子として中間流速を設定することで流れ場の変数である圧力と速度を分離する. この中間流速は上式 (3.) の圧力項をはずした式であり, 仮の流速であることに注意しなければならない. 上式 (3.) から中間流速を以下のように定義する. Δ τ Re ~ (3.3) また式 (3.) の発散を取ると Δ new new P τ Re (3.) これに式 (3.) を代入し new P Δ Δ ~ Re τ (3.5) 式 (3.) から (3.3) をひくと次式のようになる. new new P Δ ~ (3.6)
9 3--3 重み付き残差法ここで重み付き残差方程式を導く. 重み関数は, の関数であり p v, o, を用いて表す. 境界値は与えるものであり誤差はないため重み関数も零とする. 式 (.5) の場合 : S S e o e Π Π o oτ η η δ η δ (3.7) 式 (3.3) の場合 : Δ τ v v v v Re ~ (3.8) 式 (3.5) の場合 : Δ P new ~ p p (3.9) 式 (3.6) の場合 : Δ P new new ~ v v (3.) ここでガウス グリーンの定理 () Γ Γ g f fgn g f を利用して式 (3.8),(3.9) に適用すると Γ Δ Γ n v τ vτ v v v Re Re ~ (3.) Δ Γ Γ P n P new new ~ p p p (3.) このようにして導いた式を有限要素方程式に変換することで微分方程式から連立方程式へと変換することができ, 得られた連立方程式を解けばよい.
3-- 連立方程式の解法本研究では領域内を要素数 365, 接点数 653 の分割で有限要素解析を行った. 有限要素法の各処理から得られた多元連立一次方程式を繰り返し計算することで発達した流れを観察することができる. 本研究では反復法の前処理付き共役勾配法を用いて数値計算を行った.
3- プログラムの流れ開始条件設定要素分割境界値設定 τ ~ の計算の計算 本研究で作成したプログラムの流れを図 3. に示す. 従来のソフトではコーシーの運動方程式にニュートン流体の構成方程式を代入した式, すなわちナビエ ストークス方程式について数値解析を行っていたため非ニュートン流体の数値解析を行うには上で記述した重み付き残差法など, 有限要素法の手順を初めから行う必要があった. そこでこのプログラムでは偏差応力テンソルをコーシーの運動方程式と連立させて数値解析を行っている. こうすると流体が変わった場合でも, 偏差応力テンソルの部分のみを考慮すれば非ニュートン流体についても数値解析を行えるような構造になる. これがこのプログラムの特徴である. P の計算の計算 各値出力 ゼロクリアー 収束判定 終了 図 3. フローチャート
第 章計算条件と要素分割 - 計算条件 本研究では図. のような 3 次元の血管における血液の流れについて解析を行った. 図. は血管の中心部を通る - 平面の断面図である. 本研究では太い動脈を対象とした血液流動の解析を行うので, 表. より, 代表速度である流入の平均速度 5[mm/s], 代表長さである血管の直径 6[mm], 球部分の半径 3[mm], 血管の長さ 3[mm] で解析を行った. 境界条件として壁面での滑りなし, 流入部は平均速度 5[mm/s] のポアズイユ流れ, 流出部で 3 は対流流出条件, 血液の物性値より. [ P s] Re 8 として解析を行った. 3 η, ρ 5 [ kg m ], レイノルズ数 表. (3) ヒトの体循環 血管 直径 [cm] 平均速度 [cm/s] レイノルズ数 太い動脈.-.6-5 -85 毛細血管.5-..5-..7-.3 太い静脈.5-. 5- -57 D z D D 3D 図. 血管の中心を通る - 平面
3 - 要素分割 本研究では図. のような流路において要素数 365, 接点数 653 に分割して有限要素解析を行った. 座標軸の原点は流入部の円の中心にとった. z 図. 要素分割図
第 5 章結果と考察図 5. の (),(b) はそれぞれ血管の中心を通る - 平面におけるニュートン流体とキャッソン流体の速度ベクトル図である. 図 5. の (c) は (),(b) のキャッソン流体とニュートン流体の速度差のベクトル図であり, 速度差の最大値はニュートン流体とキャッソン流体をそれぞれ求めたときの最大値の約 分の となっている. 図 5. の ()~() は を 分割したそれぞれ -z 平面におけるキャッソン流体とニュートン流体の速度差のベクトル図である. これらのベクトル図からニュートン流体, キャッソン流体ともに狭窄部後方において渦が発生し逆流していることがわかるが,(c) よりキャッソン流体はニュートン流体より狭窄部によって渦が発生し逆流させられる影響が少ないと言える. すなわち狭窄部後方下部においてキャッソン流体のほうがニュートン流体より出口に向かう流速が大きいということである. これはキャッソン流体の粘度がニュートン流体の粘性を上回るためにキャッソン流体の流れがニュートン流体の流れに比べ粘性に支配されることが原因であると考えられる. () ニュートン流体 (b) キャッソン流体 (c) キャッソン流体とニュートン流体の速度差 v v 図 5. 速度ベクトル図 :- 平面 c n
5 () -. (b) -.3 (c) -. () -. (e).
6 (f). (g). ().3 (). 図 5. キャッソン流体とニュートン流体の速度差 v v :-z 平面 c n
7 次に示す図 5.3 が流線図である. 図 (),(b) は全体の流線図,(c),() は z における - 平面の流線図である. 本研究においては Re8 としたため流れは慣性にほぼ支配されている, そのため, ニュートン流体とキャッソン流体では顕著な違いが確認できなかった. new.5 -.5 3 5-.5.5 z () ニュートン流体 cs.5 -.5 3 5-.5.5 z (b) キャッソン流体
8 cs3.5 -.5 3 5 (c) ニュートン流体.5 new3 -.5 3 5 () キャッソン流体 図 5.3 流線図
9 図 5. の (),(b) はそれぞれニュートン流体, キャッソン流体の下部壁面せん断応力分布図である. 図 5. の (c) は下部壁面せん断応力の差を表している. キャッソン流体は全体的にニュートン流体より高い値をとり, どちらも狭窄部前半部分で最も高い値を示す. これはその部分は流れの剥離点の近傍であることから何らかの相関関係があるものと考えられる. 応力差については狭窄部中央付近においてもっとも差が大きくなる. 応力差の最大値はせん断応力をそれぞれ求めたときの最大値の約 5 分の となり, また, 前述したように速度差の最大値についてはニュートン流体, キャッソン流体の速度をそれぞれ求めたときの最大値の約 分の となっている. しかしこの値が大きいか小さいかの判別は対象が人間であるために最大限の注意が必要であり, それは今後の課題である.
() ニュートン流体 (b) キャッソン流体 (c) キャッソン流体とニュートン流体の応力差 τ τ c n 図 5. 下部壁面における壁面せん断応力
第 6 章結言狭窄部を有する血管内流れにおいてニュートン流体とキャッソン流体の流動解析を行い, それぞれ速度, 応力にどのような影響を与えるのかを調べ以下の結論を得た. () キャッソン流体の非ニュートン性は血管内の流速に影響を与えていた. それは狭窄部後方において著しく違いが現れた. () 壁面せん断応力分布はキャッソン流体の方が全体的に高い値をとり, それは特に狭窄部中央付近において顕著である.
参考文献 () 有限要素法による流れのシミュレーション (998) 川原睦人らシュプリンガー フェアラーク東京 () 血液のレオロジーと血流 (3) 菅原基晃, 前田信治コロナ社 (3) バイオレオロジー (98) 岡小天掌華房
3 謝辞本研究を行うにあたって, 終始懇切丁寧な御指導頂いた蝶野成臣教授ならびに辻知宏助教授に深く御礼申し上げます. また知能流体力学研究室の撰隆文氏, 田口圭一氏にも多大なご援助を頂きました. 重ねて深く御礼申し上げます.