第 24 回数値流体力学シンポジウム
B4-5

患者個別モデルに対するマルチスケール血流解析
Multi-Scale Blood Flow simulation for a Patient-Specific Model⃝

佐藤 友喜，東大院工, 東京都目黒区駒場４−６−１，E-mail：yukisato@iis.u-tokyo.ac.jp
梁 夫友，理研, 埼玉県和光市広沢２−１, E-mail：fyliang@riken.jp
高木 周, 東大工, 東京都文京区本郷７−３−１，Email：takagi@mech.t.u-tokyo.ac.jp
大島 まり, 東大情報学環, 東京都目黒区駒場４−６−１, E-mail：marie@iis.u-tokyo.ac.jp

Yuki SATO, Dept. of Mechanical Engineering, The Univ. of Tokyo, 4-6-1 Komaba Meguro-ku Tokyo
Fuyou LIANG, Integrated Simulation of Living Matter Group, RIKEN,2-1 Hirosawa, Wako, Saitama
Shu TAKAGI, Dept. of Mechanical Engineering, The Univ. of Tokyo, 7-3-1 Hongo Bunkyo-ku Tokyo
Marie OSHIMA, Interfaculty Initiative in Information Studies, The Univ. of Tokyo, 4-6-1 Komaba Meguro-ku Tokyo
In this study, we achieved to develop multi-scale blood flow simulation program that takes account of circulatory system for patientspecific geometry. Patient-Specific geometry is extracted from medical image, which is CT (computer tomography) angiography
or MRI (magnetic resonance imaging). This 3D simulation program has the boundary conditions, which are determined from 1-0dimensional(1-0D) blood flow analysis for cardiovascular system. This program controls not only boundary conditions in the outlet
but also those in the inlet of 3D model. In addition, we investigated how peripheral resistances aﬀect 3D flow state in order to know
the characteristics of peripheral resisntances. Herein we describe 3-1-0D coupling method and present the influence of multi-scale
simulation and the peripheral resistances of arteries.

1. はじめに循環器系疾患の一つである脳梗塞あるいは心筋梗塞は，動脈硬化症によって引き起こされることが多い．動脈硬化の発症メカニズムは完全には解明されてはいないが，血流による血管壁への力学的刺激が原因の一つとして考えられている．力学的刺激としては，壁面せん断応力によるずり刺激，血圧による伸展刺激や圧刺激が挙げられる
が，これらの情報の直接的な in vivo 計測には課題が多く，数値シミュレーションによる血管病変の予測や発症メカニズムの解明が期待されている．
壁面せん断応力の分布は，血管形状に起因した流れ場の影響に依存することから，数値解析では，CT や MRI 等の医用画像から血管形状を抽出し，Patient-Specific な血管モデルを用いた解析を行うことが必要とされている．また，
in vivo 環境に近い生体環境を再現するためには，血管形状に加えて，全身の血管網や毛細血管等が与える血流抵抗，血管の弾性による圧力脈波の伝播の影響を考慮することが必要とされている．特に，スケールの小さい動脈の血流抵抗を計算する場合は，医用画像の空間分解能の限界から，血管形状の抽出を行うことが困難であるため，異なる計算手法を用いて，3 次元流体解析と coupling を行う
マルチスケール解析が必要となる (1) ．Vignon-Cremental
ら (2) は，複数の分岐を持つ腹部大動脈の境界条件に電気回路を模擬した Impedance Model(0 次元モデル）を与え，末梢血管の影響を考慮した計算を行っている．また，徳田 (3) は，脳血管を解析対象として，軸方向流れのみの解析を行う 1 次元モデルと，0 次元モデルを末梢血管の血流抵抗として 3 次元解析の出口境界条件としたマルチスケール計算を行っている．しかし，現実の生体内現象としては，末梢血管の血流抵抗の変化によって，3 次元解析の血流量の流入波形の変化が起こることが想定されるため，流入波形を固定して計算を行っている点に課題が残る．そこで，全身の血液循環を考慮した計算を行い，3次元モデルの流出境界条件だけでなく，流入境界条件を算出することが望まれている．
全身の血管網と血液循環の計算を行った研究においては，Liang ら (4) (5) が，全身 55 本の動脈に対して 1 次元，細動脈や毛細血管，静脈に対して 0 次元の解析を適用し，心臓モデルを内装した全身の血液循環解析を行っている．
本研究では，3 次元血流解析に対して，1 次元-0 次元モデルによる全身の血液循環モデルを入口・出口境界条件として用いたマルチスケール解析によって，実際の生体内の血流場の再現を目標としている．

また，本研究では，末梢血管系に対して Lumped Parameter model(0 次元モデル) が用いられているが，特に 3 次元モデル下流において，これらのパラメータは血流場に大きな影響を与えることが想定される．特に，3 次元モデルより下流に存在する Lumped Parameter は，血流場に大きな影響を与えることが想定される．本研究では，全身の血液循環という境界条件の影響に加えて，末梢血管抵抗と流れ場の関連性とその傾向について考察を行う．

2. 解析手法
2.1 3-1-0 次元連成解析の概要
本研究では，左頸動脈分岐部に対して 3 次元流体解析を行い，境界条件を全身 55 本の主要な動脈に対しては 1次元解析，残りの毛細血管，静脈，心臓に対しては 0 次元モデルを用いて解析を行った．図 1 は，3-1-0 次元連成解析の概要を示す．それぞれの解析手法に関しては次節に示す．

Fig. 1: Overview of cardiovascular system
2.2 数値解析手法
2.2.1 3 次元モデル
血管壁を剛体とした流体部のみの解析を行う．支配方程式は，以下の連続の式と Navier-Stokes 方程式となる．
∇·u=0(1)∂u∇P+ (u · ∇)u = −+ ν∇2 u
∂t
ρ

(2)

u は速度ベクトル，ρ は密度，P は圧力，ν は動粘性係数，t は時間を表す．解析手法は，非構造格子を適用した有限体積法，離散化アルゴリズムは SMAC 法を用いた．
1c 2010 by JSFM
Copyright ⃝第 24 回数値流体力学シンポジウム
B4-5

2.2.2 1 次元モデル
血管壁の弾性から生じる圧力脈波の伝播の影響を考慮した軸方向の流れ解析を行う．基礎方程式には，次の形の連続の式，Navier-Stokes 方程式，圧力の状態方程式を用いる．
∂A ∂Q+=0∂t∂z( 2)∂Q∂ QQA ∂P++ 8πρµ = 0+∂t∂z Aρ ∂zA√ A
P − P0 = β − 1A0[iv] [ii],[iii] の計算を繰り返し行い，圧力 Pin(3D) の拍動
1 周期の平均値と，前回の coupling で得られた圧力の平均値の差が十分小さくなったところで計算を終了する．
3 次元-1 次元モデル間のインターフェイスにおける計算では，Lumped parameter model に準じた以下の方程式を用いる．
Qu(1D) − Qu(3D)
dPu
=(9)
dt
Cu
Qd1(3D) − Qd1(1D)
dPd1
=(10)
dt
Cd1
Qd2(3D) − Qd2(1D)
dPd2=
(11)
dt
Cd2

(3)(4)(5)

Pu − (Pin(3D) − Pout1(3D) ) − Pd1
dQu(3D)=dt
LudQd2(3D)Pd1 + (Pout2(3D) − Pout1(3D) ) − Pd2
=
dt
Ld2

ここで，P0 は参照圧力，A は血管の断面積，A0 は血管の基準断面積を表す．β は血管壁の弾性を表すパラメータを表し，次の式によって求める．

Eh0
β=r0 (1 − σ2 )(6)(12)(13)µ は粘性係数，ρ は密度， E は血管壁のヤング率，h0 は基準壁面厚さ，r0 は血管の基準半径．また，σ はポアソン比を表し，本研究では 0.5 とした．
計算手法は 2 段階 Lax-Wendroﬀ 法を用いた．
2.2.3 0 次元モデル
0 次元モデルは，主に毛細血管，静脈の解析に対して用い，図 2 で示すような Lumped parameter model を適用する (6) ．血管の弾性，血流の慣性力，粘性力をそれぞれパラメータ C,L,R で表し，式 (7),(8) を用いて，流量 Q と圧力 P の関係を求める．

Fig. 3: Algorithm
ここで，Cu ,Cd1 ,Cd2 は血管壁のコンプライアンスを表し，Lu ,Ld2 は流体のインダクタンスを表す. 式 (9),(10),(11)は，1 次元モデルと 3 次元モデルの流量変化から両モデル間の絶対圧力を算出する方程式である．ここで求めた圧力
Pu ,Pd1 ,Pd2 は，次の運動量の釣り合いを表す式 (12),(13)に代入し，3 次元モデルの流量 Qu(3D) ,Qd2(3D) を算出する．
また，3 次元モデルと 0 次元モデルの間に interface として 0 次元モデルを導入することにより，1 次元モデルの境界における圧力脈波の後進波の修正を行い，3 次元モデルの領域に存在する血管弾性の影響を考慮した 1 次元モデル末端での流量境界条件を与えることが出来る．1次元モデルでは，対流項の影響が支配的なため，粘性項を 0 と仮定すると，式 (3),(4),(5) から次の疑似線形双曲型の方程式が導出される (7) ．

Fig. 2: Lumped parameter modeldPi+ Qi+1 − Qi = 0dt(7)dQi+1= − (Pi+1 − Pi ) − RQi+1dt

(8)

CL

2.3 連成手法
本研究では，3 次元解析と 1 次元解析を拍動周期ごとに交互に計算して連成を行う．計算アルゴリズムを図 3，
3 次元-1 次元解析の coupling 部分のモデルを図 4 に示す．
計算の流れは次のようになる．

∂U∂U+ H(U)=0∂t∂x
ここで，U と H(U) は  A 
U =  Q0 √H(U) =  βA ( Q )2−2ρ A0A[i] 1 次元-0 次元モデルのみの全身循環解析を行い，得られた流量情報 Qu(1D) ,Qd2(1D) を 3 次元モデルの流入，流出口に与える．
[ii] 3 次元モデルの解析を行い，得られた圧力差の情報
Pin(3D) − Pout1(3D) ,Pout2(3D) − Pout1(3D) を 1 次元-0 次元解析に与える．
[iii] 3 次元-1 次元モデル間の Interface の計算を含めた 3次元モデル以外の残りの血管の流れ解析を行い, 算出された流量情報 Qu(3D) ,Qd2(3D) を 3 次元モデルの流入，流出口に与える．
2(14)

(15)
1 
2Q A(16)c 2010 by JSFM
Copyright ⃝第 24 回数値流体力学シンポジウム
B4-5

次元,1 次元モデル間の Lumped paramter model の計算によって求まる時刻 n における境界の圧力値を表す. これより，1 次元モデル境界の流量 Q(1D) は，An+1 ,W1 n+1 ,W2 n+1より，次のように求まる．

Q(1D) = An+1

W1 n+1 + W2 n+12(24)以上のように計算することにより，3 次元モデルの剛体壁の影響を受けずに，脈波の伝播を考慮した 1 次元モデル境界の流量を算出することが可能となる．

2.4 末梢血管モデル
総頸動脈 (CCA) から分岐する内頸動脈 (ICA) と外頸動脈 (ECA) の先には，図 (5) のように血管網を構築している．血管網の図から，Rc と R0 の抵抗値の大きさが ICAと ECA の流量分配を決める大きな要素となる．
ここで，Rc は, 式 (20) の脈波速度 c1 を用いて，次の式から求められる (9) ．
Rc =ρc1A(25)

これより，R0 の値のみが流量分配に対して支配的となることから，このパラメータを適正な値に調整する必要がある．本研究では，RT = Rc + R0 とおき，RT を変化させることで，動脈の末梢血流抵抗 R0 を変化させ，計算を行う．

Fig. 4: Coupling Interfaceを表す．ここで，H(U) を対角化すると，次に示す対角型の式が表れる．
    ∂ W1  λ1  ∂ W1 (17) +  =0∂t W2  λ2  ∂x W2 W1,2 は，特性線 λ1,2 上で不変となるリーマン不変量であり, 進行波と後進波を表す．W1,2 と λ1,2 は次のように表される．
()Q± 4 c1 − c1,0AQ
λ1,2 = ± c1
A

W1,2 =(18)(19)

Fig. 5: Peripheral Network

c1 は A = A, c1,0 は A = A0 での脈波速度を表し，それぞれ√ ( )1/4β A
c1 =(20)
2ρ A0√β(21)c1,0 =2ρ

3. 解析対象・条件
66 歳男性の左総頸動脈 (L.CCA）を解析対象とする．血管形状は 160 枚の CT 画像 (512 × 512 pixel，解像度 0.468pixel，スライス間隔 0.5 mm) より三次元再構築した．血管断面積はそれぞれ CCA：35.9 mm2 ，ECA：17.6 mm2 ，ICA：6.85 mm2 である．この血管モデルに対して，市販の計算格子作成ソフトウェア ANSYS ICEM CFD 12.0 を用いてテトラ-プリズム複合格子を作成した．要素数は254,086，節点数は 110,683 である．血液の物性値は密度を ρ = 1, 000 kg/m3 ，粘性係数を µ = 4.0 × 1.0−3 Pa · s とした．と表される．
ここで，時間 tn+1 における境界の進行波 W1 (tn+1 , L) を以下に示すような外挿によって求める (8) ．

W1 (tn+1 , L) = W1 (tn , L − λ1 ∆t)(22)

L は 1 次元モデル境界の空間座標を表す．
次に，式 (5) を変形した式(An+1 = A0)2Pn+1βFig. 6: CCA(23)

4. 解析結果

より，境界の断面積 An+1 を求め，W1 ,W2 の関係式 (18) にAn+1 から境界の後進波 W2 n+1 を求める．ここで，Pn は，334.1 マルチスケール境界条件による流量波形の変化
全身の血液循環を伴うマルチスケール境界条件を適用した解析では，流量の流出波形のみならず，流入波形も変化する．本研究では，次の 2 つの条件で計算を行い，その結果の比較を行った．
c 2010 by JSFM
Copyright ⃝第 24 回数値流体力学シンポジウム
B4-5(i) 流入境界：１次元-0 次元による全身循環解析によって算出された流量波形
流出境界：圧力 P=0 Pa
Tab. 2: Peripheral resistance 2(ii) 流入・流出境界：マルチスケール境界条件
末梢血管抵抗：(ICA : RT = 42.3 mmHg · s · ml−1 ,ECA : RT = 42.3 mmHg · s · ml−1 )RT [mmHg · s · ml−1 ] (ICA)
RT [mmHg · s · ml−1 ] (ECA)

Default

Case4

Case5

Case6

42.3
42.3

42.3
63.2

42.3
84.042.3104.5

図 (7),(8) は，それぞれの結果の流量波形を示す．
14CCAECA
ICA

12Flow rate [ml/s]

108
14Qu(3D) (CCA)
Qd1(3D) (ECA)
Qd2(3D) (ICA)6124
10
Flow rate [ml/s]

202

2.2

2.4

2.6

2.8

3Times [s]Fig. 7: Flow rate waveform (i) (B.C. p=0)

86420

1415

Qu(3D) (CCA)
Qd1(3D) (ECA)
Qd2(3D) (ICA)

1215.215.415.6
Times in 3D calculation [s]

15.8

16−1

Fig. 9: Flow rate (Case1:RT =34.0 mmHg · s · ml )Flow rate [ml/s]

10864201515.2

15.4

15.6

15.8

1614Times [s]Fig. 8: Flow rate waveform (ii) (Multi-Scale B.C.)

10Flow rate [ml/s]

この結果から，全身の血管網や血液循環の影響により，流入波形が変化していることが分かる．しかし，(ii) のマルチスケール計算では，時間 t = 15.4 s 付近で，ICAと ECA の流量の大きさが逆転する現象が見られた．この現象は，ICA 側末梢側の血流抵抗が大きすぎる，またはECA 末梢側の血流抵抗が小さすぎるために起きた可能性が考えられる．次節より，末梢血管抵抗が流量に与える影響について説明する．

015Case2

Case3

34.0
42.3

29.8
42.3

25.7
42.3

15.2

15.4
15.6
Times in 3D calculation [s]15.8

16−1Fig. 10: Flow rate (Case2:RT =29.8 mmHg · s · ml )14Qu(3D) (CCA)
Qd1(3D) (ECA)
Qd2(3D) (ICA)

1210Flow rate [ml/s]

Case1

42.3
42.3

62

Tab. 1: Peripheral resistance 1Default

844.2 末梢血管抵抗が流量・流れ場に与える影響
ICA 末梢の抵抗を表 (1), (2) のように変化させたとこ
ろ，図 (9)〜(14) のような流量波形が得られた．Case １〜3 では，ICA 側の末梢血管抵抗 RT を小さくするとともに，基準値 RT = 42.3 mmHg · s · ml−1 で見られた ICA 側とECA 側の流量の大きさが逆転する現象は見られなくなった．また，Case4〜6 では，ECA 側の末梢血管抵抗 RT を大きくするとともに，流量の逆転する現象は見られなくなった．これより，ICA 側の末梢の血流抵抗を ECA 側の抵抗に対して，相対的に小さくする必要があることが明らかになった．

RT [mmHg · s · ml−1 ] (ICA)
RT [mmHg · s · ml−1 ] (ECA)

Qu(3D) (CCA)
Qd1(3D) (ECA)
Qd2(3D) (ICA)1286
4201515.2

15.4
15.6
Times in 3D calculation [s]

15.8

16−1

Fig. 11: Flow rate (Case3:RT =25.7 mmHg · s · ml )4c 2010 by JSFM
Copyright ⃝第 24 回数値流体力学シンポジウム
B4-5
1410

Qu(3D) (CCA)
Qd1(3D) (ECA)
Qd2(3D) (ICA)

12Qu(3D) (CCA)
Qd1(3D) (ECA)
Qd2(3D) (ICA)

8Flow rate [ml/s]

Flow rate [ml/s]

10866442201515.2

15.4
15.6
Times in 3D calculation [s]15.80165−1

Fig. 12: Flow rate (Case4:RT =63.2 mmHg · s · ml )14202530
RT (ICA) [mmHg.s.ml-1]

35104045Qu(3D) (CCA)
Qd1(3D) (ECA)
Qd2(3D) (ICA)

8108

Flow rate [ml/s]

Flow rate [ml/s]

15Fig. 15: Averaged flow rate (RT (ECA) = 42.3 mmHg·s·ml−1 )Qu(3D) (CCA)
Qd1(3D) (ECA)Qd2(3D) (ICA)

121064264201515.2

15.4
15.6
Times in 3D calculation [s]15.816040−1Fig. 13: Flow rate (Case5:RT =84.0 mmHg · s · ml )
1460

7080
RT (ECA) [mmHg.s.ml-1]

Flow rate [ml/s]

RT [mmHg · s · ml−1 ] (ICA)
RT [mmHg · s · ml−1 ] (ECA)
Flow rate [mL/s] (CCA)
Flow rate [mL/s] (ECA)
Flow rate [mL/s] (ICA)
Flow rate division [%] (ECA)
Flow rate division [%] (ICA)

8642015.2

15.4
15.6Times in 3D calculation [s]15.8100110Tab. 3: Comparison between default and modified cases

101590Fig. 16: Averaged flow rate (RT (ICA) = 42.3 mmHg·s·ml−1 )Qu(3D) (CCA)
Qd1(3D) (ECA)
Qd2(3D) (ICA)

1250

16Fig. 14: Flow rate (Case6:RT =104.5 mmHg · s · ml−1 )また，図 (15),(16) は，それぞれ ICA,ECA の末梢血管
抵抗 RT を基準値 RT = 42.3 mmHg · s · ml−1 から変化させた場合の流量・流量分配の変化を表す．図 (15) では，CCA，ICA の流量のみに変化が現れた．また，図 (16) も同様に，CCA,ECA の流量のみに変化が現れた．これより，ICA，ECA の絶対流量は，ICA,ECA それぞれに対応した末梢血管抵抗に依存することがわかる．
実際の生体の流量・流量分配の関しては，近年，MRIによる血管内の流量・流速データの計測が行われている．
Marshall らは，実際の生体内の流量・流量分配を，MRIによって 14 本の頸動脈内の流量・流量分配の計測を行っている (10) . その結果によると，各血管での 1 拍動周期における平均流量は，それぞれ CCA： 6.16 ml/s, ECA：1.59 ml/s, ICA： 4.14 ml/s であった．この計測データと本研究で行った計算結果を比較すると，ICA ： RT = 16.7mmHg · s · ml−1 , ECA ： RT = 42.3 mmHg · s · ml−1 の場合の計算が計測データに最も近い値となる．
末梢血管の血流抵抗 RT を上記の値に設定し，基準値の値で設定した場合での血管内の壁面せん断応力，流れ場等の結果の比較を行った．表 (3) は比較対象の計算条件と流量を表す．

Default

Modified

42.3
42.3
3.47
1.56
1.86
45.7
54.5

16.7
42.3
5.76
1.55
4.21
26.973.1

ここで，血管壁にかかる 1 拍動周期当りの壁面せん断応力の平均値 TAWSS(Time-Avraged Wall Shear Stress) を比較したところ，図 (17) のような結果が得られた．TAWSSは次式で表される，∫1 T
T AWS S =(26)|t s | dt
T 0

ECA 側の流量は変化がなく，ICA 側の流量が増加したために，ICA 側の血管壁の TAWSS が増加したものと思われる．
ま た ，拍 動 に よ る 壁 面 せ ん 断 応 力 の 影 響 を ，OSI(Oscillatory Shear Index) に よって 確 認 し た と ころ，図 (18) のような結果が得られた．OSI (OscillatoryShear Index) とは，血流の拍動に伴い血流が血管壁に及ぼす壁面せん断の方向が空間的に変化する割合を示す指標であり，次式で表される (11) (12) ．
1(τmean )OS I =
1−(27)
2
T AWS S5

c 2010 by JSFM
Copyright ⃝第 24 回数値流体力学シンポジウム
B4-5

ここで，τmean は１拍動 T [s] における壁面上のトラクションベクトル t s の時間平均値であり，次の式で表される．
∫1 Tt s dt(28)τmean =T 0

OSI は，動脈硬化症との関連が強く指摘されている指標である．この結果より，ICA 側の流量の増加により，分岐直後の ICA 壁面近傍での壁面せん断応力の振動が少ないことが明らかとなった．
また，両ケースの流れ場を，Default の場合においてICA と ECA 側の流量の逆転が起きた t = 15.45 s 付近で可視化したところ，図 (19) のような結果が得られた．
Default の場合では，分岐直後の ICA 側に再循環流が確認された．この現象が，ICA 側の血管壁に OSI の大きな分布を発生させた一因であると考えられる．

(a)

(b)
Fig. 18: OSI(a)Default: RT (ECA)= 42.3 mmHg·s·ml−1 ,RT (ICA) = 42.3 mmHg · s · ml−1 ). (b)Modified: RT (ECA)= 42.3 mmHg · s · ml−1 , RT (ICA) = 16.7 mmHg · s · ml−1 ).

(a)

(b)

(a)

Fig. 17: TAWSS(a)Default: RT (ECA)= 42.3 mmHg ·−1s · ml , RT (ICA) = 42.3 mmHg · s · ml−1 ). (b)Modified: RT(ECA) = 42.3 mmHg · s · ml−1 , RT (ICA) = 16.7 mmHg · s ·ml−1 ).

(b)Fig. 19: Stream line (a)Default: RT (ECA)= 42.3 mmHg ·s · ml−1 , RT (ICA) = 42.3 mmHg · s · ml−1 ). (b)Modified: RT(ECA) = 42.3 mmHg · s · ml−1 , RT (ICA) = 16.7 mmHg · s ·ml−1 )., t = 15.45 s

6c 2010 by JSFM
Copyright ⃝第 24 回数値流体力学シンポジウム
B4-5

5. おわりに本研究では，全身の血液循環を考慮した左頸動脈分岐部の血流解析を行い，以下の結論を得た．
(1) 血液の全身循環により，流出境界条件だけではなく，流入境界条件に対しても流量波形の周期的な変化が起こることを確認した．
(2) ICA，ECA に流れる流量は，それぞれの血管の末梢血管抵抗に依存することが明らかになった．
(3) ICA の末梢側の血流抵抗を減少させた場合，CCA，ICA の流量が増大し，ICA の壁面近傍に発生する再循環流の減少に伴う，OSI の減少を確認した．
今後は様々な血管モデルを用いてシミュレーションを行
い，超音波ドップラー計測や MRI 計測によって得られる
in vivo の流速データを用いて，末梢血管抵抗が与える流量への影響を調査し，抵抗と流量の相関関係をさらに明らかにする．また，超音波ドップラー計測や MRI 計測によって得られる in vivo の流速データとシミュレーションの結果をより詳細に比較し，流れ場を実際の生体環境内に近づけるような末梢血管の Lumped parameter の傾向について考察する．
謝辞
本研究は，文部科学省最先端・高性能汎用コンピュータの開発利用「次世代生命体統合シミュレーションソフトウェアの研究開発」の支援を受け行われたものである．
ここに記して関係者各位に謝意を表する．
参考文献(1) C. A. Taylor, et al., Ann. Biomed. Eng., 38,3,1188-1203(2010) .

(2) I. E. Vignon-Clemental et al., Comput. Methods Appl.
Mech. Engrg. 195,3776-3796 (2006) .
(3) 徳田茂史ら, 日本機械学会論文集（B 編), 74 巻 740
号 (2008-4).
(4) F. Y. Liang, et al., J. Biomech.,42,692-704 (2009).
(5) F. Y. Liang, et al., Med. Biol. Eng. Comput.,47,743-755(2009).
(6) V. Milisic, et al., Math. Model. and Numer. Anal.,
38,613-632 (2004) .
(7) M. A. Fernandez, et al., Multiscale Model. Simul.
4,1,215-236 (2005)(8) L. Formaggia, et al., Comput. Meth. Biomech. Biomed.
Eng., 9,5,273-288 (2006)(9) G. Atlas, et al., Cardiovasc. Eng., 9,11-17 (2009)(10) I.Marshall et al., Physiol. Meas., 25,691-697(11) D. N. Ku, et al., Arteriosclerosis, 5,293-302 (1985)(12) C. A. Taylor, et al., Ann. Biomed. Eng., 26,975-987(1998)

7c 2010 by JSFM
Copyright ⃝