第 26 回数値流体力学シンポジウム
講演番号 D01-2

流体-構造連成手法を用いた患者個体別頸動脈の解析
Analysis of atient-specific carotid artery bifurcation using fluid-structure interaction method○ 片桐賢吾, 芝浦工大院工, 東京都江東区豊洲３－７－５, E-mail：md11027@sic.shibaura-it.ac.jp
石上雄太, 東大院工, 東京都目黒区駒場４－６－１, E-mail：gami@iis.u-tokyo.ac.jp
大島まり, 東大情報学環, 東京都目黒区駒場４－６－１, E-mail：marie@iis.u-tokyo.ac.jp
山本創太, 芝浦工大, 東京都江東区豊洲３－７－５, Email：sota@sic.shibaura-it.ac.jp
Kengo KATAGIRI, Dept. of Mechanical Engineering, Shibaura Inst. of Tech., 3-7-5 Toyosu, Koto-ku, TokyoYuta ISHIGAMI, Dept. of Mechanical Engineering, The Univ. of Tokyo, 4-6-1 Komaba, Meguro-ku, TokyoMarie OSHIMA, Interfaculty Initiative in Information Studies, The Univ. of Tokyo, 4-6-1 Komaba, Meguro-ku, Tokyo
Sota YAMAMOTO, Dept. Engineering Science & Mechanics, Shibaura Inst. of Tech., 3-7-5 Toyosu Koto-ku Tokyo
The hemodynamic factors such as wall shear stress are known to play an important role in vascular diseases, andrecent studies have indicated that the arterial-wall deformation have an influence on them. We carry out numericalfluid structure interaction analyses to investigate the influences of material property of arterial wall on blood flow
behavior and wall shear stress. The simulation is performed for the carotid artery, in which both lumen & arterial wallis constructed from CT images. The results show a significant difference in stress/strain and wall shear stress
distributions.

し，流体領域の低壁面せん断応力の分布にどのような影響を与えるのかを検証する．

1．緒言
重篤な循環器系疾患の原因となる動脈硬化症等の血管病変に，壁面せん断応力 (Wall Shear Stress：WSS)，伸展刺激，圧刺激といった血管壁への力学的刺激が関与していると考えられている．特に，低WSS説(1)が広く知られていて，有力視されている．これらの血行力学的影響と病変のメカニズムを解明するために，CTやMRIといった医用画像診断装置から得られた画像データにより患者独自の血管形状をモデリング化し，さらに数値シミュレーションを行う image-based modeling and simulationと呼ばれる手法が取り入れられている．
この手法はTaylorら(3), (4) が頸動脈や腹部大動脈の解析に用いており，バイオメカニクスの分野で主流になりつつある(2)．また，実際の診断や治療と深く関わりながら開発されている．従来は，血管内腔内の血流（流体）解析，あるいは血管壁の構造解析のように独立した数値解析による研究が多かったが，近年は計算機性能の飛躍的な向上により，血流と血管壁の相互作用を考慮した解析が可能となってきている(5)- (7)．X. L. Yangら(8)は肺動脈を，鳥井ら(9) は脳動脈を対象とした解析を行っている．V.
Chabannesら(10)は2次元と3次元の解析に適応した高次の流体‐構造連成解析アプリケーションの開発を行っており，様々な径の血管に大変形を伴う連成解析を適用した結果を発表している．invivo simulationを行う際には，より現実に近い生体環境を模擬する必要があり，血管形状のみならず，血管壁の材料モデルや血流，あるいは血管壁に対する境界条件の数理モデル化が重要な課題である．
本研究では動脈硬化症の発生・進行メカニズムを解明することを目的とし，血管壁の材料特性，あるいは境界条件の違いが解析結果に与える影響を検証し，考察する．WSS は血管内の流動構造に依存するため，血管の弾性変形や，流出境界条件が流量分配に与える影響を正確に把握する必要がある．手法としては，医用画像から血管内腔と血管壁の形状をモデリングし，生体環境を数理モデル化した境界条件および材料モデルを用いた in vivo 流体構造連成（Fluid-Structure Interaction）解析を有限要素法にて行う．材料モデルは健常な頸動脈，病変した頸動脈，脳血管の 3 種の非線形材料特性を与え，血管壁の弾性変形の挙動やひずみ分布を比較

2． 解析手法
2．1 流体領域における支配方程式
解析対象である頸動脈の管径は血液中の赤血球の大きさに対して十分に大きく，非ニュートン性の影響は無視できるため，血液は非圧縮性 Newton 流体とみなす．流体領域における支配方程式
は連続の式(1)と Navier-Stokes 方程式(2)である．

u i
0xi

f(1)

 ij
u iu  f cj i t
x j
x jc j  u j  u mj

(2)

(3)

 u i ij   p ij    x ju j 
xi (4)ここで，u は流体の速度，um は計算格子の移動速度，c は相対速度，ρf は流体の密度，σ は流体の表面力，p は流体圧力，μ は粘性係数，δ はクロネッカーのデルタである．
2．2 構造領域における支配方程式
血管壁は応力とひずみの関係が非線形な特性を有する．本研究では構造領域に超弾性体材料の材料特性を与えている．構造領域の支配方程式は平衡方程式である．

1

0s v  S F T (5)S ij (6)W Eij EijCopyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 D01-2
ここで，v は構造の変位ベクトル，ρs は構造の密度，F 変形勾配テンソル，S は第2 Piola-Kirchhoff 応力テンソル，E はGreen-Lagrange ひずみテンソルである．

W  C1 I C  3  C 2 II C  3  C3 I C  32

 C 4 I C  3( II C  3)  C5 ( IIc  3) 2  C 6 ( I C  3) 3

 C 7 ( I C  3) 2 ( II C  3)  C8 ( I C  3)( II C  3) 2
 C9 ( II C  3) 3  C10 ( I C  3) 4  C11 ( I C  3) 5(7)
ここで，I1，I2 は右Cauchy-Green 変形テンソルC の第1，第2 主不変量である．非線形材料特性は(7)式のひずみエネルギー密度関数
W により与えられ，C1 ～ C11 の定数により実験値とフィッティングさせる．
2．3 離散化手法および連成手法
流体の支配方程式にALE （Arbitrary Lagrangian-Eulerian) 有限要素法，構造の支配方程式に混合型有限要素法を用いて離散化する．
連成では図1に示されるような強連成に基づく．

Fig. 2 Schematic view of multi-scale boundary pressure analysisq   4 q 2  A pq   8t x  3 A   x
Ap ( x, t )  p 0 

時間積分法に関しては，構造領域ではNewmark-β法を，流体領域ではNewmark-β法を拡張したPMA（Predictor Multi-correctorAlgorithm）法を使用している．流体解析と構造解析の連成には，Galerkin法により式(1)，(2)，(5)から得られた有限要素方程式を重ね合わせ，境界面上の幾何学的連続条件(8)および平衡条件(9)を満たすよう強連成系を構築して，Newton-Raphson法により解く．

u it ij n j  Tij n j  0Z  x,   (12)P  x,  Q  x,  (13)

ここで，Z はインピーダンス，ω は角速度である．
2．3 解析モデル
解析対象は動脈硬化症の好発部位とされる左総頸動脈(Lt.
Common Carotid artery：CCA)，内頸動脈(Lt. Internal Carotid artery：ICA)および外頸動脈(Lt. External Carotid artery：ECA)の分岐部周辺とする．血管形状は66歳男性のCT画像(512×512pixel, 解像度0.468mm/pixel, スライス間隔0.5mm) 160枚から3次元形状へ再構築したものである．血管形状を図3に示す．血管径は2.9 ～ 12.3 mm，壁厚はCCA で0.5 ～ 1.2 mm，ICA で0.5 ～ 1.6 mm，ECA で0.8～ 1.0 mmであり，血管長は74.9 mmである．計算格子は6面体1次要素を生成し，解析モデル全体で51,720節点，流体領域で 32,090要素，構造領域で 16,696 要素を有し，構造領域は半径方向へ4層を生成している．

(8)(9)

ここで，n は流体と構造の境界面の法線ベクトル，T は構造のCauchy 応力テンソルである．
2．4 末梢血管抵抗を考慮した流出境界圧力解析
本研究では流出境界に末梢血管抵抗の影響を考慮した1D－0D解析(11)から得た圧力波形を与えた．3D－1D－0Dモデルのカップリングによるマルチスケール解析の概要を図2に示す．
1Dモデルは解剖学的知識に基づくフラクタル性ネットワーク構造(12) , (13)を有する，3次元モデルの下流から血管径が100 μmに至るまで対称に分岐する動脈，小動脈にあたる血管を想定している．
支配方程式は連続の式(10)，Navier-Stokes方程式(11)および圧力と断面積の関係式(12)である．

q A
0x tA0 
A 

ここで，A は断面積，Q は流量，P は圧力， E はヤング率，h0 は基準壁面厚さ，r0 は初期血管半径，σ はポアソン比である．
細動脈，毛細血管にあたる，血管径が100 ～ 5 μmまでをSTI（Structure Tree Impedance）モデル(14)を用いた0Dモデルにて解析し，1D解析モデルでの流出境界条件としてカップリングする． STIモデルのインピーダンスは，線形化された1D解析の支配方程式に基づき，波数空間にて表現されている．したがって，圧力は流量と関連する方程式(13)にて表される．

Fig. 1 Schematic view of strong coupled methodvi 4 Eh 1
3 r0 (11)3． 解析条件
3．1 物性値
流体領域の作動流体はヒトの血液の物性値を想定しており，密度は1,000 kg/m3, 粘性係数は4.06 x 10-3 Pa・secである．構造は生体材料を模擬するため，超弾性体材料の高次Mooney-Rivlin体を想定し，(7)式における11個の材料定数を表1のように設定することで健常な総頸動脈（Model A）(15)，粥状プラークを有す総頸動脈（ModelB）(15)，脳血管（Model C）(16)の3種の材料特性を図3の血管壁領域に与える．図4に使用した材料特性の応力－ひずみ曲線を示す．

(10)2Copyright © 2012 by JSFM第 26 回数値流体力学シンポジウム
講演番号 D01-2
4． 解析結果および考察
本研究にて使用する1D-0D解析による圧力境界条件の影響を，自由流出境界条件と比較する．材料Cの解析に2つの境界条件を与えた時の流出境界ICA，ECAでの圧力波形を図6に示す．1D-0D解析と連成することで，ICAとECAの間の圧力差が減少し，波形が図5の流入波形に近付いていることがわかる．また，波形の位相遅れも自由流出条件に比べ小さくなっている．このとき，流入流量を100とした場合の流量分配比をICA : ECAで表すと，1D-0D境界条件は75 : 25，流出境界条件は80 : 20であり，圧力境界条件の差異が流量分配に影響を与えていることがわかる．
次に，1D-0D境界条件を与えた解析に3つの材料モデルを適用した場合のICA，ECAにおける圧力波形を図7に示す．ICA，ECAともに材料C，B，Aの順に平均的に圧力が高い．これは想定する材料特性が固くなるほど血圧が高くなる結果といえる．また，流入波形のピークは62 ms，ICAの流出波形のピークは91 msに，ECAの流出波形のピークは88 msに現れており，材料モデルによる差異はみられない．このとき，流入流量を100とした場合の流量分配比
は，材料Aで74: 26，材料Bで74 : 26，材料Cで75 : 25であり，顕著な差異は現れなかった．これらの結果より，流出境界の圧力波形と流量分配は1D-0D解析の結果に大きく依存し，材料モデルによる差異は小さいものと考えられる．

ICA

CCA

Pressure

Flow

ECA

Fig. 3 Patient-specific carotid artery bifurcation geometry
Tab. 1 Mooney-Rivｌin constants
C value

1234567891011A

6.7105.220043.900014.4

1.33B

2.1603.16002320002.78

3.07C

17501000010010000000020←peak at inlet

Case C ICA (1D-0D)
Case C ECA (1D-0D)
Case C ICA (free-stream)

Pressure [ x 10 3 Pa]

Case C ECA (free-stream)

16

Fig. 4 Stress-Stretch curves of material model A(healthy common
carotid artery), B(diseased common carotid artery) and C(cerebralvessel artery)[15],[16]12これら応力－ひずみ曲線に基づき，ヒトの一般的な血圧である80～ 120 mmHgの範囲で適用している．
3．2 境界条件および初期条件
流体領域の流入境界には拍動流入条件を適用し，流量の時間に対するプロファイルは，20代男性の総頸動脈からエコー計測された拍動波形を17次関数でフィッティングし，Poiseuille分布として与えている．計測データと算出した拍動波形を図5に示す．流出境界条件は1D-0Dモデルにより算出された時間に対する圧力プロファイルとして与えている．ここで，ICA下流では 1Dモデル17分岐，0Dモデル14分岐を有し，ECA下流では1Dモデル16分岐，0Dモデ
ル14分岐を有している． 構造の境界条件は流入･流出端面が移動・変形しないものと仮定し，境界端面の節点変位を完全拘束している．また，初期条件は流体部の流速と圧力，および構造部での応力をゼロとする．

80.00.20.40.6
Time [s]0.81.0Fig. 6 Pressure wave form at ICA and ECA in Case C.

19

Case A ICA
Case A ECA
Case B ICA
Case B ECA
Case C ICA
Case C ECAPressure [ x 10 3 Pa]←peak at inlet

15110.00.20.40.6
Time [s]0.81.0Fig. 7 Pressure wave form at ICA and ECA. Case A is healthy
common carotid artery material, Case B is diseased common carotid
artery material, Case C is cerebral blood vessel material.

Fig. 5 Time-flow rate wave form (measurement and approximated)
3Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 D01-2
1拍動中で流速が最大となる時刻（72 ms）における，各材料モデルの構造領域における最大主ひずみの分布を図8に示す．この結果は1拍動中でひずみの大きい領域の運否する傾向はすべての材料モデルで共通しているものの，材料A，Bと比較して材料Cのひずみは1/10程度であり，非常に固い材料特性であるといえる．材
料A，B，Cの解析結果におけるひずみの最大値はそれぞれ0.038，0.034，0.0027である．また，構造領域の最大変位の数値も材料A，B，Cでそれぞれ0.89 mm，0.76 mm，0.073 mmであり，これらから脳血管の材料特性は総頸動脈の材料特性と比較して顕著に固い挙動を示すことが明らかである．
最大主ひずみ分布と同時刻の，構造領域における最大引張主応力の分布を図9に示す．材料B，Cでは応力値の高い領域が大きく，特に材料Cでは分岐近傍で非常高い応力域が見られる．しかし，材料AではICA側の分岐部近傍で局所的に高い応力を示し，Bの最大応力値を超えている．これは生成した計算格子による影響も考えられるため，今後より詳細な検証を行う必要がある．

Model A
Max = 0.038

WSS[ Pa ]

Model C
Model B
Model A
Fig. 10 Wall shear stress distribution
1 拍動中で流速が最大となる時刻（72 ms）における壁面せん断応力の分布を図 10 に示す．
動脈硬化症は低壁面せん断応力領域で生じる傾向にあることが知られ，工藤ら(17)の実験ではプラークの原因となる LDL（Low Density Lipo-protein）の血管壁への取り込みは 0 ～ 4 Pa の領域において増大するとしている．よってここ
では拍動ピークの結果に対し WSS = 0 ～ 4 Pa の範囲に着目し可視化を行っている．すべてのモデル分岐部の外側に低 WSS の分布が見られ，材料 A では血管全体で低 WSS せん断応力域が広範囲に分布しており，材料特性が固くなるほど低 WSS 領域が縮小される傾向にある．またどの材料モデルにおいても分岐部周辺で低 WSS 域が発生している．特に材料 C の結果では ICA の分岐直後に流れの停滞による低 WSS 域が局所的に発生しており，これはこの患者の右総頸動脈分岐部の動脈硬化症による血管狭窄部位と一致している．これにより，血管壁の材料特性が固くなることで血管全体の低 WSS 域は縮小されるものの，より局所的な低WSS の発現と，それによる病変の発生が起こりうることが示唆される．

Model BMax = 0.034

Max principle
Strain [ - ]Model C

5． 結言
本研究では，CT画像から再構築した総頸動脈分岐部の血管内腔および血管壁形状モデルに対して種々の非線形材料特性および境界条件を適用してFSI解析を行い，以下の結論を得た．
(1) 1D-0D解析モデルの導入により流出境界間の圧力波形の差が減少した．また材料特性の差異によりECA境界で圧力波形に差異が生じたが，流量分配に顕著な影響を与えなかった．
(2) 構造領域のひずみ，応力分布は同一の傾向を有する者の，材料特性の差異による影響がひずみ，応力，変位において顕著に確認された．
(3) 血管壁の材料特性が固いほど低WSSの領域は減少したが，同時に病変好発部位に局所的な低WSS領域が発現した．

Max = 0.0027

Fig. 8 Max principle strain distribution

Model A Max = 638 Pa

謝辞
本研究で使用した医用画像は藤田保健衛生大学の早川基治先生よりご提供いただいた．ここに記して関係者各位に謝意を表する．

Model B Max = 548 Pa

Max principle
Stress [ Pa ]Model C

参考文献
(1) Malek. A. M., et al., Hemodynamic shear stress and its role inatherosclerosis. JAMA, Vol. 282, pp. 2035–2042, (1999)(2) C. A. Taylor, et.al., Image-Based Modeling of Blood Flow and
Vessel Wall Dynamics: Applications, Methods and Future
Directions, Annals of Biomedical Engineering, vol. 38, no. 3, pp.
1188-1203, (2010)(3) C. A. Taylor, et al., Finite element modeling of blood flow in arteries,
Computer Methods in Applied Mechanics and Engineering, Vol 158,

Max = 1125 Pa

Fig. 9 Max principle stress distribution
4Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 D01-2
pp. 155- 196, (1998)(4) C. A. Taylor, et al., Patient-Specific Modeling of Cardiovascular
Mechanics, Annual Review of Biomedical Engineering, Vol. 11:
109-134, (2009)(5) K. Takizawa, et al., Fluid–structure interaction modeling andperformance analysis of the Orion spacecraft parachutes,
International Journal for Numerical Method in Fluids, Vol. 65,( 1-3):pp.271-285 (2011)(6) T. J. R. Hughes, et al., Isogeometric fluid–structure interaction
analysis with applications to arterial blood flow, Computational
Mechanics Journal, Vol. 38:pp.310–322, (2006)(7) Y. Bazilevs, et al., 3D simulation of wind turbine rotors at full scale.
Part II: Fluid–structure interaction modeling with composite blades,
International Journal for Numerical Method in Fluids, Vol. 65,( 1-3):pp.236-253, (2011)(8) X. L. Yang, et al., Fluid-structure interaction in a pulmonary arterial
bifurcation, Journal of Biomechanics, Vol 40, (12): pp. 2694-2699,(2007)(9) Ryo Torii, et al., Influencing Factors in Image-Based Fluid-StructureInteraction Computation of Cerebral Aneurysms, International
Journal for Numerical Method in Fluids, Vol. 65, ( 1-3):pp.324-340,(2011)(10) V Chabannes, et al., High order fluid structure interaction in 2D and3D application to blood flow in arteries, Journal of Computationaland Applied Mathematics, ISSN 0377-0427, (2012)(11) A Krdey, et al., Fluid-Structure Interaction Simulation of Middle
Cerebral Artery Using Multi-Scale Model as Outflow condition, 第23 回バイオエンジニアリング講演会講演論文集, pp.533-534,(2011)(12) S. Tokuda, et al., Modeling of Outflow Boundary Conditions
Taking Account of the Effects of Peripheral Vessel Network for
Hemodynamic Simulation, Transactions of the Japan Society ofMechanical Engineers, B, 74 (740): pp. 871-878, (2008)(13) Zamir, M., FractalDimensions and Multifractility in VascularBranching, Journal of Theoretical Biology, Vol. 212, (2): pp.
183-190, (2001)(14) Olufsen, M. S., et al., Numerical Simulation and Experimental
Validation of Blood Flow in Arteries with Structured-Tree Outflow
Conditions, Annals of Biomedical Engineering, Vol. 28, 1281-1299,(2000)(15) H. Yamada, et al., Correlation between Stress/Strain and theRetention of Lipoproteins and Rupture in Atherematous Plaque ofthe Human Carotid Artery, Journal of Biomechanical Science andEngineering, Vol. 5( 4): pp. 291-302, (2010)(16) 林紘三郎, バイオメカニクス, コロナ社, (2000)(17) S. Kudo, et al., Effect of Wall Shear Stress on Macromolecule
Uptake into Cultured Endothelial Cells, Transactions of the Japan
Society of Mechanical Engineers, B, 64 (618): pp. 367-374, (1998)

5Copyright © 2012 by JSFM