— 3 —SENAC Vol. 49, No. 1（2016. 1）［共同研究成果］

拍動流下における動脈瘤を有する血管壁挙動の数値解析法とその臨床応用
西川優、横堀壽光、工藤和貴、大見敏仁、門間良平、友野雄基、菅野崇史(東北大学大学院工学研究科)

1. 緒言
拍動流によって誘起される血管壁の動的挙動の解明には流体解析と流体圧力により誘起される血管壁の変形挙動の構造解析を行うことが必要である．
流体解析には，2 次元あるいは軸対称流れにおいては渦度―流れ関数法，3 次元流れに対しては MAC 法がある．前者は，連続の式を厳密に満足するが圧力を求める解析への適用が困難である．[1] 一方，MAC 法は連続の式が近似的にしか満足されていないことが指摘されている．[1]ところで，著者の一人は以前より超音波ドップラー効果を利用して，拍動下での血管粘弾性発現度と拍動の不規則度を定量化して動脈硬化と動脈瘤を非侵襲で診断する手法を提案してきた．[2-5]特に，動脈瘤の存在予測診断法の確立には，血管壁の不規則拍動挙動の測定実験結果を数値流体力学（CFD）により検証することが重要である．そこで，著者らは，渦度―流れ関数法と圧力のポアソン方程式を連成して解く方法により，動脈瘤を有する血管壁の拍動挙動を数値解析により明らかにしてきた．[6-8]
本稿では，これらの解析手法について，さらに著者のグループで検討・発展させてきた手法とそれにより得られた結果について，2 節で使用する動脈瘤モデルと解析の流れ，3節で基礎式と差分方法，4 節で具体的な解析，5 節で結言の順に述べる．

2. 動脈瘤モデル
本研究で用いた解析モデルを Fig.1 に示す．動脈瘤形状は(1) 式に示す Gauss の正規分布曲線を用いて表した．[11]R0 は初期管半径，b,c はそれぞれ動脈瘤高さ，幅を決定するパラメータである．

� � �(z) � �� � ��� exp(−�� ��� �)(1)

また，流体は均質な非圧縮性 Newton 流体，流れは軸対称流とし，流れ方向および半径方向の二次元で解析を行う．
次に，解析領域を壁面形状関数で無次元化することにより，(R, Z)座標系から(η, Z)座標系(長方形領域)への変換による写像を施す． (1) 式で表される動脈瘤壁面形状関数における独立変数を無次元化すると， (2) 式で表される[11]．

— 4 —

SENAC Vol. 49, No. 1（2016. 1）

�(�) � � � �����(���� )(2)(2) 式による(R, Z)座標系から(η, Z)座標系への変換は式(3)で与える．

��������� � �
�(�)

(3)

これを図 2 に示す.

Fig.1 解析モデル

1v

u0
Fig.2

(η, Z)座標

2.1 解析のながれ
本解析のフローチャートを図 3 に示す．
はじめにΨ，Ω の初期値を設定しΨについて SOR 法で解く．それを用いてΩについて陽解法で解く．このΩを用いてまたΨを求め以下それを繰り返す．
次に圧力 P については Navier-Stokes の方程式から圧力のポアソン方程式を求める．

P の増分と現在の血管径 rn 及び有限変形の式[9]から計算される剛性を用いて rn＋１を求める．
また安定的に計算を行うため，一定ステップ数定常流を流し，十分落ち着いた時点で拍動流に切り替える．

— 5 —

拍動流下における動脈瘤を有する血管壁挙動の数値解析法とその臨床応用

Fig.3 フローチャート

3. 基礎式
二次元円筒座標系における連続の式と Navier-Stokes 方程式を (4) ,(5) 式に示す．

∂�� 1 �(�� � ) = 0+
� ��∂z(4)1 ��� � �� 1 ��� � � �����
�������
�
�+�+
�
=−+�
�++�� �� �� ��� �����z�� ��� ��

�1 ��� � � � 1 �� � � � � � � ��� ��� ����+ ��=−+�� � ++− ��� � + ��� ��� �����z���� ��� ��u+ :

軸方向速度[m/s]v+ :半径方向速度[m/s]

+t :

時間[s]

p:

圧力

[Pa=N/m2]ρ:

密度[kg/m3]

ν:

動粘性係数

[m2/s](5a)(5b)

— 6 —

SENAC Vol. 49, No. 1（2016. 1）

添え字の+は有次元量であることを表す．ここで (4) 式より式(6)を満たす Ψ が存在する．

�� = −� ��

�� = −

� ��� ��(6)� ��

また，軸対称流において渦度ベクトル ω は (7) 式で定義される．

�=�� � ���−����(7)

次に一様管部の半径を R0，平均流速を U0 で表し無次元数 R, Z, Ψ, Ω,u,v,t,P を導入する．
�=

���

，� =

���

，Ψ =���� ����，Ω =u=��� �����
��，v =

����

，t =

�� ��
� ，P =
���� �����

ここで U0 は(8)式のレイノルズ数の定義から，Re を定めることにより決定した．

�� =

�� ��
�

(8)

以上の無次元変数および無次元流れ関数 Ψ，無次元渦度 Ω を導入することにより(4)，(5)式から以下の Poisson 方程式および渦度輸送方程式を導出される．
(Poisson 方程式)

�� =(渦度輸送方程式)
����+���
��+����������� −�− Ω=�� ��� ������+(9)����� �������+�� ��� �+��� ��(10)�� �

また，境界条件は以下のように定める．

管入口:

管出口:

管中心軸 :

壁面上:���� ��������

��= 0,= 0,�� ����

����= 0,�� ����

=0= 0, P = 100������

= 0, v = 0, Ψ = 0, Ω = 0,

u = 0, v = 0, Ω =定常流：Ψ = −1,� �� � ��,� ��� ��

��

��

=0

振動流：Ψ = −1 × (2 − cost),

=0

(11a)
(11b)
(11c)

(11d)

— 7 —

拍動流下における動脈瘤を有する血管壁挙動の数値解析法とその臨床応用

壁面上の流れ関数に関する境界条件は定常時は-1 の一定値とし，拍動時には全域で cos 関数に従って変動し振動を再現している．
なお，一階微分，二階微分は (12)-(15) 式で表される.

1 �∂=∂R �(�) ��(12)1∂�∂�=∂R� (�(�))� ∂��(13)� ��(�) ∂∂∂=−+�(�) �� ∂� ∂�
∂Z(14)

������(�) ∂�
2η��(�)� � � �(�) ∂∂�=�
�
+��
�
−��
��
�(�) �� � ∂�∂Z � ��(�)�� ��∂��
��(�)�� ��(�) ∂�∂�
−2+�(�) �� ∂�� �� ∂� �(15)

これらを基礎方程式(9)(10)に適用すると，無次元基礎方程式は以下のように表される．
(Poisson 方程式)

(渦度輸送方程式)
����=�������� ������� ����� + ��+ ���� ���� +����
��+ ���� ��� ���� ��� ��+ ������+ ���� ��� �(16)= �� Ω+ ��� Ω� − �������− �������− ��� Ω(17)各係数を式(18)に示す．
�1
��(�)
�� =
�1 + �� �� ����(�(�))
�11
��(�)� � �(�)
�� =
�−+2�
�
�
−��(�)����� �(�(�))�
�
�� = −2� ��(�)�(�) ��(18a)
(18b)
(18c)

�� = 1(18d)�� = 1(18g)�� = ��(�)�� = ��

�� = ��(18e)
(18f)
(18h)

— 8 —

SENAC Vol. 49, No. 1（2016. 1）

�11
��(�)
� � �(�)
�� =
�+2�
�
�
−
��(�)
����� �(�(�))� �
��� = − ���� =��� =(18i)1 ����(�)(18j)2���
�(�(�)) ����� = −(18k)2���
�(�(�)) ��(18l)2� ��(�) �� ���−��
�(�(�)) �(�) �� �� ��(18m)(16)，(17)式を差分化して解き，Ψ と Ω を計算することによって Navier-Stokes 方程式の解を求める．
得られる Ψ から，軸方向速度 u と半径方向速度 v は式(19)のように求めることができる．

�=−

�

��

�(�(�))� ��， �=���(�)�−�� ��(�) ��

�(�) ��

��+����

�

(19)

流れ関数-渦度法では圧力を求める事ができないため，ここまでの計算より得られるΨ，u,v を用いて圧力のポアソン方程式を解き P を計算する．
(5a)式を Z で偏微分し，(5b)式に r+をかけ r+で偏微分し 1/r+をかける．それらを足しあわせる．ここで連続の式の項 D を式(20)のように定義する．

� ���� =∂�� � � �� �+ �+ �∂z�
��

(20)

連続の式を満たせば 0 となるはずであるが，差分化の際の離散化誤差が集積し大きな値となる場合があるため式(21)の形で D 項を式に残す．[1]��1 ��∂� D1 �����∂� D
�� = − � � + �� � + � � � � + � �++��������∂� � � � � �� � ∂� � �������� �� ��� ���− �( � )� + 2 � � + ( � )� + � ����� ����
��(21)

— 9 —

拍動流下における動脈瘤を有する血管壁挙動の数値解析法とその臨床応用

3.1 差分化Poisson 方程式に二次精度中心差分，渦度輸送方程式に時間微分に一次精度上流差分，空間微分に二次精度中心差分を用いて離散化する．
(Poisson 方程式)
�������
���Ψ������� − (��)� Ψ����� −

�

�

��

��
� �2 �(��)� + (��)� � Ψ��� + (��)� +�������

�����
��

���

Ψ������� = �� Ω�����Ψ������� + �(��)� −� Ψ����� −������������� Ψ����� −��Ψ������� + (��)� Ψ����� +(22)(渦度輸送方程式)�Ω������ − Ω���

Δt=1����Ω�− 2Ω���� + Ω������ ��� (Δ�)� �����+���Ω�− 2Ω���� + Ω������ �(Δ�)� �����+����� ��Ω������ − Ω������ � + ��� Ω���� � −�Ω− Ω������ �2Δ�2Δ� �����+−���Ω�− Ω�������� − Ω�������� + Ω�������� ��Δ�Δ� �������(23)����Ω�− Ω������ � − ��� Ω����
2Δ� �����(22) 式に対して SOR 法を適用し流れ関数Ψを導出する．そこで得られたΨの値を用いて(23)式により Euler 陽解法を用いて渦度 Ω を求める．なお (21) 式を SOR 法にて計算する際の収束条件は以下のように定め，加速緩和係数 β は 1.5 とした．
(���)�� ���(�)�����(�)����(���)����� � 1���(�)

(m: 計算ステップ数)(���)�= ���� + � ��
���

(�)

− ���� �

� : Gauss-Sidel 法による解，β:加速緩和係数)(�
流体解析と同様に圧力式についても写像と差分化を式(24)に示すように行う．

— 10 —

���� =SENAC Vol. 49, No. 1（2016. 1）

0.5S(Z)� ℎ�

1 + S(Z)� + ����(�)������� ������ + ��������
��(�)
ℎ�

������ − ������2� � ��� −� �� ��
�(�)
�(�)
2ℎ+�−2� � �������� − �������� − �������� + ���������
�(�)4ℎ�+1 ������ + �������(�)�
ℎ�+������ + ������1 ������ − ������+���(�)�
ℎ2ℎ+ �� �−(24)����� ��� ���++ � �−��+��(�) �� ��Δt�(�) ��2� � �������
�� 2� � � � � ��
�� ��
�+�−
�−
�−�
�(�)�
�(�)
�
�(�)� �������
��+���
1 ��1 ���++��� � ��(�)� �� �(�)� ���

12� � �� ��2 �� ���� ���+ �� �−�+
�(�)� ���(�)� �� �� �(�) �� ��+ �−����� �� ���+ � + ����(�) �� ��� �(�)�

ここで�� =である．

��(�)��� �� =� � �(�)�� �

圧力についても SOR 法を用い，収束条件は下の式とした．
(���)�� ���(�)�����(�)����(���)����� � 10��(�)

(m: 計算ステップ数)(���)= ���� + �������(�)− ���� �(��: Gauss-Sidel 法による解，β:加速緩和係数)

— 11 —

拍動流下における動脈瘤を有する血管壁挙動の数値解析法とその臨床応用

3.2 圧力による血管の変形
血管壁が管内圧力に応じて変形する効果を考える．内圧 p における血管内径は薄肉円管モデルにおける Laplace の式[11]を用いて， 漸化式とし(25) 式で表される．
������ = ������
��
��

(25)

その様子を図 4 に示す．ここで各パラメータは以下のように定義した．
E:

管壁の剛性

[Pa=N/m2]h:

管壁の厚さ[m]p

Pn -Pn-1[Pa]rn 今計算ステップの管半径rn 1 前計算ステップの管半径[m][m]Fig.4 各ステップの圧力と半径の関係
ここで管内圧力の変動により得られた血管形状で次ステップの流れ場の計算を行い，そこから再び圧力の計算および血管形状の計算を行う．このように数値流体力学による流れ場の計算と材料力学による血管壁形状の計算の連成解析を行うことで血管内での現象を再現している．管内圧力に変動が生じることに伴い， (25) 式より r に変動が生じるが，このr に対して逐次式(3)で表される(R, Z)座標系から(η, Z)座標系(長方形領域)への変換による写像を施すことにより，拍動流条件下における血管の脈動を一様流条件下と同一の手法で解析することが可能となる．
ここで， (25) 式に示す管壁の剛性は，Fung による理論式(26)より[9]σ=�∗��∗ � ��∗��� � � )������� � �∗ ))

(26)

�σ:

血管壁に作用する応力λ:ひずみ[g/mm2]から式(27)を計算することにより算出した．

E=����

=�∗��∗ � ��∗���� � � � �� � � �� ������� � �∗ ))
��

(27)

ここで a，σ*，λ*は定数であり，体重 15～25kg のイヌの胸部大動脈を用いた引張強度試験による実験結果から[10]，以下の値を用いた．
a:

1.9

*

σ :

25.83

λ* :

1.6

[g/mm2]

— 12 —

SENAC Vol. 49, No. 1（2016. 1）

4. 解析
著者の一人による過去の研究により，動脈瘤を持つ人の血管壁速度波形には，図 5(b)のように頂点にピークが 2 回表れる二相性と呼ばれる特徴が表れることが知られている[4]．測定により二相性の有無を判定することで動脈瘤を診断する手法を提案しており[4]，CFD によりその検証が行われている．[6-8]そこで本解析では二相性が動脈瘤由来であるかを原理的に確かめるため，上記の解析手法を用いて計算した血管壁速度波形に臨床と同様の二相性が表れるのかを検証した．また二相性の特徴を表す，図 6 のように波形の時間発展を一つの図で表すことのできる血管壁拍動波形の軌跡解析が行われている．1 ループが 1 拍動であり，軌道は長周期的なずれを生じる．健常者は計測定例では周期的なループを描き，それに対して動脈瘤を有する場合は，
図 7 に示すように 2 相性に対応する小さな circuit loop が右上部に発生する[4]．
本論文の CFD 解析により動脈瘤を有する血管壁の振動挙動に臨床波形測定例と同様の二相性波形が発現するか検証する．

(a) 健常者(b)動脈瘤患者
Fig.5 血管壁速度波形例[4](a) 血管壁速度測定結果(b)アトラクターによる拍動軌跡
Fig.6 波形軌道解析例[6,12]

X(t+Δt)

X(t+Δt)拍動流下における動脈瘤を有する血管壁挙動の数値解析法とその臨床応用

X(t)(a) 健常者の波計測定例— 13 —

CircuitX(t)(b)動脈瘤を有する波形測定例Fig.7 解析結果のアトラクター解析[4]4.1 CFD による解析結果
本解析では(1)式の動脈瘤形状を決定する係数 b=1.5，c=1 として解析を行った．速度波形
測定場所は図 8 に示すモデル上流部(a)と，中央部(b)である．図 9 に解析により得られた血管壁速度波形を示す．上流部(a)，中央部(b)共に動脈瘤を有する血管に臨床測定例と同様の二相性波形を再現することができた．また波形軌道解析結果を図 10 に示す．右端に circuitloop が発現している．
以上のことから拍動による血管壁挙動の周期変化において，測定による二相性波形は動脈瘤の存在により発現することが示された．また二相性は動脈瘤部から離れた上流部でも発現した．

(a)

(b)

Fig.8 壁面速度測定位置

— 14 —

SENAC Vol. 49, No. 1（2016. 1）

0.0000004

Wall Steps

Wall Speed

0.0000002
0-2E-07-4E-07-6E-07-8E-07224500

288250Time Steps

0.000006
0.000004
0.000002
0-0.000002-0.000004-0.000006-0.000008-0.00001-0.000012
224500(a)288250Time Steps(b)

Fig.9 血管壁速度波形

Fig.10 解析結果

5. 結言
本解析により非侵襲測定で得られた図 4,図 6 の二相性波形と Circuit loop は動脈瘤の存在に寄るものであることが CFD により検証された．今後実構造に近い血管壁への解析に本プ
ログラムを発展させることにより，さらに詳細に動脈瘤の存在による血管壁挙動の乱れを再現することができると考えられる．
謝辞
本解析手法の開発にあたり，東北大学サイバーサイエンスセンターの SX-ACE を使用させて頂いたことを感謝する．
参考文献[1] 河村哲也 流体解析１，1996,朝倉書店[2] A.T.Yokobori,Jr., T.Ohkuma, S.Sasak, H.Yoshinari, T.Yokobori, H.Ohuchi andS.Mori,

Bio-Medical Materials and Engineering, 4, (1994), pp.87−96.

[3] A.T.Yokobori,Jr., M. Ichiki, H.Ohuchi, T.Kobayashi, T.Satoh and Y.Kinoshita,
Bio-Medical Materials and Engineering, 14,3 (2006), pp.241−249.
[4] A.T. Yokobori,Jr., M. Owa, M. Ichiki, T. Satoh, Y. Ohtomo, Y. Satoh, S. Ohogoshi, Y.

— 15 —

拍動流下における動脈瘤を有する血管壁挙動の数値解析法とその臨床応用

Kinoshita and S. Karino, Journal of Atherosclerosis and Thrombosis, 13,4,(2006),pp163-174.

[5] A.T.Yokobori,Jr.,T.Ohmi,

R.Monma,

Y.Tomono,K.Inoue,M.Owa,M.Ichiki,

N.Mochiduki and H.Yamashita, Bio-Medical Materials and Engineering, 23, (2013),
pp.75−91.
[6] A.T.Yokobori.Jr., T.Ohmi, Y.Tomono, R.Monma, Y.Nishikawa, T.Kanno M.Ichiki The
analysis of pulsatile behavior of blood vessel wall with aneurysm based on hybrid
analysis of computer fluid dynamics (CFD) and solid mechanics, Proc of annual
meeting of Atherosclerosis and Thrombosis, (2013)[7] A.T.Yokobori.Jr., T.Kanno, Y.Nishikawa, T.Ohmi, The analyses of mechanicalbehaviors of blood vessel wall with aneurysm at the progressive stage based on
simulate experiment and computer fluid dynamics(CFD), Proc of annual meeting ofAtherosclerosis and Thrombosis, (2014)[8] A.T.Yokobori.Jr., Y.Nishikawa, K.Kudo, T.Ohmi, Long Periodic Fluctuation of Blood
Vessel Wall Caused by the Viscoelasticity Under Pulsatile pressure Flow, Proc ofannual meeting ofAtherosclerosis and Thrombosis, (2015)[9] Y.C.Fung, Appl. Mech. Rev., 1968,21(1),pp1-20.
[10] 横堀壽光，大熊恒郎，前山俊秀，横堀武夫，大内博，奈良茂樹，佐々木久雄，葛西森夫，生体材料，1986,4,pp27-32.
[11] 岡小天，レオロジー入門，工業調査会，1970.
[12] D. Ruelle, and F. Takens, Communications in Mathematical Physics, 1971, 20,
pp167-192.