Microsoft Word 卒業論文2.doc

平成 16 年度

卒業論文

狭窄部を有する血管内の血流の有限要素解析

高知工科大学

工学部

知能機械システム工学科

知能流体力学研究室

清水 昌彦

１－２ 血液の性質 1 １－３ 数値計算 1 １－４ 有限要素法の概要 2 第２章 基礎方程式 ２－１ 支配方程式 4 ２－１－１ 連続の式 5 ２－１－２ コーシーの運動方程式 6 ２－１－３ 血液の構成方程式 6 ２－２ 無次元化 7 第３章 解析手法 ３－１ 有限要素解析 8 ３－２ 手順 8 ３－２－１ 時間方向の離散化 8 ３－２－２ 流速修正法 8 ３－２－３ 重み付き残差法 9 ３－２－４ 連立方程式の解法 10 ３－３ プログラムの流れ 11 第４章 計算条件と要素分割 ４－１ 計算条件 12 ４－２ 要素分割 12 第５章 結果と考察 14 第６章 結言 21 参考文献 22 謝辞 23

第１章 緒言

１－１ 本研究を行う背景と目的 現在，死亡原因の中の多くには心疾患や脳血管疾患がある．これらの病気の発生メカニ ズムは血液流動と血管に深く関わっており，疾患と血液流動との関係を解明することが望 まれている．しかし，実際のヒトの血管を用いて実験することは困難であり，流体力学の 豊富な知識が必要とされるため，医学分野のみで血液流動を解析するのは不可能である． よって，流体力学の立場から考えることが必要である．そこで，近年バイオメカにクスの 分野では人体実験をする必要のないコンピュータシミュレーションによる解析は今まで行 われてきたが，血液をニュートン流体として解析されたものがほとんどで，血液の非ニュ ートン性を考慮したものは少ない．そこで本研究では血液の非ニュートン性を考慮して血 液流動の解析を行う．解析対象として半球のついた円管内流れとし，血管狭窄の状態に近 づけた．この状況下でのニュートン流体，非ニュートン流体の流動を解析し比較する． １－２ 血液の性質 血液は酸素や二酸化炭素の運搬や，栄養物や老廃物の運搬などヒトにとって重要な役割 を果たす． ヒトの体内を循環する血液は約 45%の細胞成分と約 55%の液体成分から構成される．細 胞成分には赤血球，白血球，血小板があり，液体成分は血漿である．血液に含まれる赤血 球の容積比率は約 45%と非常に大きい．よって血液流動には赤血球の変形が大きく関わっ ている．血液は水のようなニュートン流体ではなく，血液に与えるせん断速度を大きくす れば見かけ粘度が小さくなる非ニュートン流体である．さらに血液の通り道である血管は 粘弾性管である． 血液の非ニュートン性を表す式にキャッソンの式がある．キャッソンは印刷インクの流 動挙動についてずり速度とずり応力の関係を円錐－平板型粘度計で測定して式を導いたが， この式が血液流動にも使えることがわかった． １－３ 数値計算 数値計算とは解析的に解くことが不可能な自然現象などをモデル化し数値的に近似解を求 め，自然現象などを可視化することができる方法である．対象によって数値計算の方法も 使い分けなければならない．本研究においては血液流動の解析を行うため，実際の血管形 状や，疾患部を考慮し有限要素法を使うことにする． 流動の解析を行う際，境界で変化が激しいことから有限要素法は計算領域に非構造格子 を用いるため境界付近で格子を小さくし，精度よく解析することが可能である．

１－４ 有限要素法の概要 ある微分方程式 − f =0 dx du f = const( .>0) in 0≤ x≤1 (1.1) がある．近似解を u とすると f dx u d r= − (1.2) となり， r は残差である．これにある関数 * u をかけて積分したものを 0 にしなければなら ない．よって，

∫

Ω Ω 0= * rd u (1.3) となる．関数_φ*_{は重み関数といい，重み関数によって得られる方程式を重み付き残差方程} 式という．有限要素法では各要素を独立したものとして扱うので各要素によって近似関数 を定義しなければならない．各要素の近似関数は接点値を元に直線で近似する．この近似 関数を要素補間関数という．図 1.1 のように長さ h の要素に着目し，接点値を満たすように 要素補間関数を定義すると b a u h x u h x u ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ − = 1 (1.4) ここで ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ − h x 1 ， ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ hx は形状関数という． a

u

u

b

u

a b h 図 1.1 要素補間関数

また，重み関数についても同様に近似関数を導入すると * * * 1 _a u_b h x u h x u ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ − = (1.5) となり，式(1.4)，(1.5)を式(1.3)に代入すると ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ − + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛− ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ −

_∫

∫

1 0 1 0 * 1 1 1 1 a b a dxu h h x dxu h h x u 0 1 1 1 0 1 0 * ₌ ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛− ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + b

∫

a

∫

dxub h h x dxu h h x u (1.6) 式(1.6)を積分すると 0 2 1 2 1 2 1 1 2 1 1 2 2 * 2 2 * = ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛− + ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ − + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛_{− +} b a b b a a u h u h u u h h u h h u (1.7) ここで重み関数をu*a =1， 0 * ₌ b u とすると 0 2 1 1 2 1 1 2 2 ⎟_⎠ = ⎞ ⎜ ⎝ ⎛ − + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛_{− +} = _{a b} a u h h u h h N (1.8) また

u

_a*

=

0

，

u

_b*

=

1

とすると 0 2 1 2 1 2 2 ⎟_⎠ = ⎞ ⎜ ⎝ ⎛ + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛− = _{a b} b u h u h N (1.9) ここでN ，_a N はそれぞれ_b u ，*_a u によって求められたことを表すために付け加えた．b* 式(1.8)， (1.9)をまとめると ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ − − + − = ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ b a b a u u h h h h h h N N 2 2 2 2 2 1 2 1 2 1 1 2 1 1 (1.10)

この式(1.10)を有限要素方程式といい，この方程式は着目した要素内において接点の未知 関数に関する量

u

_a，

u

_bで離散化した方程式系となっている．この式を領域内すべての要素 において求め，有限要素方程式を領域全体で重ね合わせることで領域の支配方程式を単純 な連立方程式に変換し解く．このように微分方程式を要素ごとに近似方程式にし，領域全 体で重ね合わせ，支配方程式を近似した連立方程式を解くというのが有限要素法の原理で ある．

第２章 基礎方程式

ここではまず主な記号について説明する． ij e ：変形速度テンソル( i ， j )成分 ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ = i j j i ij x u x u 2 1 e L ：代表長さ U ：代表速度 P ：圧力 y S ：降伏応力 t ：時間 η ：キャッソン粘度 γ& ：ずり速度

( )

2e_ij ρ ：密度 ji τ ：偏差応力テンソル Π ： 2 ij e ２－１ 支配方程式 本研究において用いられている支配方程式は，連続の式，コーシーの運動方程式および キャッソン流体の構成方程式であり，式の表記法として総和規約を用いる．また，添字 i， j は整数で 1～3 の範囲を持つ．i=1,2,3のとき x，y，z 方向の方程式を表す． ２－１－１ 連続の式 連続の式は質量保存の式で流体が非圧縮のとき次式で表される． 0 = ∂ ∂ i i x u (2.1) ここで，x は座標成分であり，_i x₁=x,x₂ =y,x₃ =zである． u は速度の_i x 成分である． _i

２－１－２ コーシーの運動方程式 コーシーの運動方程式は運動量保存の式であり，領域に外部からなされる力がない場合， 次式のようになる． j ji i j i j i x x P x u u t u ∂ ∂ + ∂ ∂ − = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ τ ρ (2.2) ２－１－３ 血液の構成方程式

(

ij − Sy

)

= τ η γ& 1 (2.3) 上式(2.3)は血液の構成方程式としてよく用いられるキャッソンの式である．この式を τ に ついて整理すると次式のようになる． ij e Π Π τ 2 4 2 4 4 ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ + + = y yη η ij S S (2.4) しかし，Π=0のとき，τ が微分不可能になるという特性があるので修正キャッソンの式 を導入しΠ=0のときでも微分可能にするためにδを加えた． ij y y ij S S e Π Π τ 2 4 2 4 4 ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ + + + + = η δ η δ (2.5) δはずり速度が低い領域での血液粘度にあうように定めた．また，血液の構成方程式の降 伏応力がないとき，つまりSy =0のときはニュートン流体の構成方程式である．

２－２ 無次元化 ２－１の支配方程式を無次元化する． * i i Uu u = (2.6) * i i Lx x = (2.7) * y y S L U S =η (2.8) * 2 P U P=ρ (2.9) * t U L t= (2.10) ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ = _* * * * i j j i ij x u x u L U e (2.11) * Π Π 2 2 L U = (2.12) ここで代表長さ L には血管の直径，代表速度 U には流入の平均流速を用いた．*は無次元量 を表している． 上式(2.6)から(2.12)を式(2.5)に代入すると * * 4 * * * * * * 4 2 2 * 2 2 * 2 1 4 2 4 2 1 4 2 4 ij ij y y i j j i y y ij L U e L U S S x u x u L U L U S L U L U S L U τ Π Π Π Π τ * * * * η η δ δ η δ δ = ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ + + + + = ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ + + + + = (2.13) 次にコーシーの運動方程式(2.2)に式(2.6)，(2.7)，(2.9)，(2.10)および(2.13)を代入すると * * * * 2 * * * * * j ji i j i j i Lx L U Lx P U Lx Uu Uu t U L Uu ∂ ∂ + ∂ ∂ − = ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ _ρ η τ ρ 上式を両辺 L U2 ρ で割る，すなわち， 2 U L ρ をかけると * * * * * * * * * * * * * Re 1 j ji i j ji i j i j i x x P x UL x P x u u t u ∂ ∂ + ∂ ∂ − = ∂ ∂ + ∂ ∂ − = ∂ ∂ + ∂ ∂ τ τ ρ η (2.14) となる．(Re：レイノルズ数)

第３章 解析手法

３－１ 有限要素解析 本研究では重み付き残差法によって数値解析を行う．重み付き残差法とは微分方程式の 残差と重み関数の内積を零に近づけていく方法である．ここではその実際の手順を述べる． ３－２ 手順 ３－２－１ 時間方向の離散化 コーシーの運動方程式に対して時間方向の離散化 j ji i new j i j i new i x x P x u u t u u ∂ ∂ + ∂ ∂ − = ∂ ∂ + Δ + τ Re 1 (3.1) 0 = ∂ ∂ i new i x u (3.2) ３－２－２ 流速修正法 流速修正法とは流速の予測子として中間流速を設定することで流れ場の変数である圧力 と速度を分離する．この中間流速は上式(3.1)の圧力項をはずした式であり，仮の流速である ことに注意しなければならない．上式(3.1)から中間流速を以下のように定義する． ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ − ∂ ∂ Δ − = j ji j i j i i x x u u t u u τ Re 1 ~ _(3.3) また式(3.1)の発散を取ると ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ − ∂ ∂ ∂ ∂ − ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ − Δ = ∂ ∂ j ji j i j i i i i new i i new x x u u x x u x u t x P τ Re 1 1 2 2 (3.4) これに式(3.2)を代入し i i j ji j i j i i i i new x u t x x u u x x u t x P ∂ ∂ Δ = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ − ∂ ∂ ∂ ∂ − ∂ ∂ Δ = ∂ ∂ ~ 1 Re 1 1 2 2 _τ (3.5) 式(3.1)から(3.3)をひくと次式のようになる． i new i new i x P t u u ∂ ∂ Δ − = ~ (3.6)

３－２－３ 重み付き残差法 ここで重み付き残差方程式を導く．重み関数は x, の関数でありy o,v,pを用いて表す． 境界値は与えるものであり誤差はないため重み関数も零とする． 式(2.5)の場合： Ω + Ω ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ + + + + = Ω

∫

∫

∫

Ω Ω d Ω d S S d ij ij y y ij e oe Π Π o oτ η η δ η δ 4 2 2 2 4 4 (3.7) 式(3.3)の場合： ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ Ω ∂ ∂ − Ω ∂ ∂ Δ − Ω = Ω

∫

∫

∫

∫

Ω Ω Ω x d Ω x d u u t d u d u j ji j i j i i τ v v v v Re 1 ~ _(3.8) 式(3.5)の場合： Ω ∂ ∂ Δ = Ω ∂ ∂

∫

∫

Ω Ω _x d u t d x P i i i new ~ 1 2 2 p p (3.9) 式(3.6)の場合：

∫

∫

∫

Ω Ω Ω ∂ Ω ∂ Δ − Ω = Ω d x P t d u d u i new i new i v~ v (3.10) ここでガウス・グリーンの定理(1)

∫

Ω

∫

Γ

∫

Ω ∂ Ω ∂ − Γ = Ω ∂ ∂ d x g f d fgn gd x f i i i を利用して式(3.8)，(3.9)に適用すると ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ Ω ∂ ∂ + Γ − Ω ∂ ∂ Δ − Ω = Ω

∫

∫

∫

∫

∫

Ω Ω Ω x d Γ n d Ω x d u u t d u d u j ji j ji j i j i i v τ vτ v v v Re 1 Re 1 ~ _(3.11) Ω ∂ ∂ Δ = Ω ∂ ∂ ∂ ∂ − Γ ∂ ∂

∫

∫

∫

Γ Ω Ω x d u t d x P x d n x P i i i new i i i new ~ 1 _p p p (3.12) このようにして導いた式を有限要素方程式に変換することで微分方程式から連立方程式へ と変換することができ，得られた連立方程式を解けばよい．

３－２－４ 連立方程式の解法

本研究では領域内を要素数 346051，接点数 64531 の分割で有限要素解析を行った．有限 要素法の各処理から得られた多元連立一次方程式を繰り返し計算することで発達した流れ を観察することができる．本研究では反復法の前処理付き共役勾配法を用いて数値計算を 行った．

３－４ プログラムの流れ

条件設定

開始

各値出力

要素分割

ゼロクリアー

収束判定

終了

の計算

τ

u

の計算

P

の計算

u~

の計算

境界値設定

条件設定

開始

各値出力

要素分割

ゼロクリアー

収束判定

終了

の計算

τ

u

の計算

P

の計算

u~

の計算

境界値設定

図 3.1 フローチャート 本研究で作成したプログラムの流れを図 3.1 に示す．従来のソフトではコーシーの運 動方程式にニュートン流体の構成方程式を 代入した式，すなわちナビエ・ストークス 方程式について数値解析を行っていたため 非ニュートン流体の数値解析を行うには上 で記述した重み付き残差法など，有限要素 法の手順を初めから行う必要があった．そ こでこのプログラムでは偏差応力テンソル をコーシーの運動方程式と連立させて数値 解析を行っている．こうすると流体が変わ った場合でも，偏差応力テンソルの部分の みを考慮すれば非ニュートン流体について も数値解析を行えるような構造になる．こ れがこのプログラムの特徴である．

第４章 計算条件と要素分割

４－１ 計算条件 本研究では図 4.2 のような 3 次元の血管における血液の流れについて解析を行った． 図 4.1 は血管の中心部を通る x-y 平面の断面図である．本研究では太い動脈を対象とした血 液流動の解析を行うので，表 4.1 より，代表速度である流入の平均速度 500[mm/s]，代表長 さである血管の直径 6[mm]，球部分の半径 3[mm]，血管の長さ 30[mm]で解析を行った．境 界条件として壁面での滑りなし，流入部は平均速度 500[mm/s]のポアズイユ流れ，流出部で は対流流出条件，血液の物性値よりη=4.0×10−3

[

Pa⋅s

]

，

[

3

]

10 1050× kg m = ρ ，レイノルズ数 800 Re= として解析を行った． 表 4.1(3) ヒトの体循環 血管 直径[cm] 平均速度[cm/s] レイノルズ数 太い動脈 0.2-0.6 20-50 110-850 毛細血管 0.0005-0.001 0.05-0.1 0.0007-0.003 太い静脈 0.5-1.0 15-20 210-570 図 4.1 血管の中心を通る x-y 平面

D 3D

D

x

z

y

D

U

４－２ 要素分割 本研究では図 4.2 のような流路において要素数 346051，接点数 64531 に分割して有限要 素解析を行った．座標軸の原点は流入部の円の中心にとった．

x

z

y

_x

z

y

_x

z

y

図 4.2 要素分割図

第５章 結果と考察

図 5.1 の(a)，(b)はそれぞれ血管の中心を通る x-y 平面におけるニュートン流体とキャッソ ン流体の速度ベクトル図である．図 5.1 の(c)は(a)，(b)のキャッソン流体とニュートン流体 の速度差のベクトル図であり，速度差の最大値はニュートン流体とキャッソン流体をそれ ぞれ求めたときの最大値の約 1000 分の 1 となっている．図 5.2 の(a)~(i)は y を 10 分割した それぞれ x-z 平面におけるキャッソン流体とニュートン流体の速度差のベクトル図である． これらのベクトル図からニュートン流体，キャッソン流体ともに狭窄部後方において渦が 発生し逆流していることがわかるが，(c)よりキャッソン流体はニュートン流体より狭窄部 によって渦が発生し逆流させられる影響が少ないと言える．すなわち狭窄部後方下部にお いてキャッソン流体のほうがニュートン流体より出口に向かう流速が大きいということで ある．これはキャッソン流体の粘度がニュートン流体の粘性を上回るためにキャッソン流 体の流れがニュートン流体の流れに比べ粘性に支配されることが原因であると考えられる． (a)ニュートン流体 (b)キャッソン流体 (c)キャッソン流体とニュートン流体の速度差vc−vn 図 5.1 速度ベクトル図：x-y 平面

(a) y=-0.4

(b) y=-0.3

(c) y=-0.2

(d) y=-0.1

(f) y=0.1

(g) y=0.2

(h) y=0.3

(i)y=0.4

次に示す図 5.3 が流線図である．図(a)，(b)は全体の流線図，(c)，(d)は z=0 における x-y 平面の流線図である．本研究においては Re=800 としたため流れは慣性にほぼ支配されてい る，そのため，ニュートン流体とキャッソン流体では顕著な違いが確認できなかった． (a)ニュートン流体 (b)キャッソン流体

0

1

2

3

4

5-0.5

0

0.5

-0.5

0

0.5

x

z

y

new21

0

1

2

3

4

5-0.5

0

0.5

-0.5

0

0.5

x

z

y

cas21

(c)ニュートン流体 (d)キャッソン流体 図 5.3 流線図

0

0

1

2

3

4

5

-0.5

0

0.5

x

y

cas3

0

0

1

2

3

4

5

-0.5

0

0.5

x

y

new3

図 5.4 の(a)，(b)はそれぞれニュートン流体，キャッソン流体の下部壁面せん断応力分布 図である．図 5.4 の(c)は下部壁面せん断応力の差を表している．キャッソン流体は全体的に ニュートン流体より高い値をとり，どちらも狭窄部前半部分で最も高い値を示す．これは その部分は流れの剥離点の近傍であることから何らかの相関関係があるものと考えられる． 応力差については狭窄部中央付近においてもっとも差が大きくなる．応力差の最大値はせ ん断応力をそれぞれ求めたときの最大値の約 50 分の 1 となり，また，前述したように速度 差の最大値についてはニュートン流体，キャッソン流体の速度をそれぞれ求めたときの最 大値の約 1000 分の 1 となっている．しかしこの値が大きいか小さいかの判別は対象が人間 であるために最大限の注意が必要であり，それは今後の課題である．

(a)ニュートン流体

(b)キャッソン流体

(c)キャッソン流体とニュートン流体の応力差τc−τn

第６章 結言

狭窄部を有する血管内流れにおいてニュートン流体とキャッソン流体の流動解析を行い， それぞれ速度，応力にどのような影響を与えるのかを調べ以下の結論を得た． (1)キャッソン流体の非ニュートン性は血管内の流速に影響を与えていた．それは狭窄部後 方において著しく違いが現れた． (2)壁面せん断応力分布はキャッソン流体の方が全体的に高い値をとり，それは特に狭窄部 中央付近において顕著である．

参考文献 (1)有限要素法による流れのシミュレーション(1998) 川原睦人ら シュプリンガー・フェアラーク東京 (2)血液のレオロジーと血流(2003) 菅原基晃，前田信治 コロナ社 (3)バイオレオロジー(1984) 岡小天 掌華房

謝辞

本研究を行うにあたって，終始懇切丁寧な御指導頂いた蝶野成臣教授ならびに辻知宏助 教授に深く御礼申し上げます．また知能流体力学研究室の撰隆文氏，田口圭一氏にも多大 なご援助を頂きました．重ねて深く御礼申し上げます．