— 24 —SENAC Vol. 48, No. 4（2015. 10）[研究成果]［共同研究成果］
CFD による血管壁拍動挙動における Vasomotion の再現解析
工藤和貴、横堀壽光、西川１．

優、大見敏仁、門間良平、井上恭介(東北大学大学院工学研究科)

緒言

拍動流によって誘起される血管壁の動的挙動は、血管疾患に伴う壁の力学的性質変化により変化すると考えられ、その解明には流体解析と流体圧力により生じる血管壁の変形挙動の構造解析が必要である。
流体解析には、2 次元あるいは軸対称流れにおいては、渦度―流れ関数法、3 次元流れに対しては、MAC 法が一般的には用いられる。しかしながら、前者は連続の式を厳密に満足するが、血管壁の動的挙動の解明に必要な圧力を求める解析への適用が困難である[1]。また、MAC 法は、連続の式が近似的にしか満足されていないことが指摘されている[1]。
このような問題を解決するため、著者らは、渦度―流れ関数法と圧力のポアソン方程式を連成して解く方法により、血管壁の拍動挙動を解析している[6-8]。血管疾患に影響する血管拍動を構成する振動成分の特定は、より高精度に血管疾患を診断するために必要である。この手法により、拍動下での血管壁粘弾性発現度を定量化して動脈硬化の進行度を非侵襲で診断することが可能となる。
本稿では、著者のグループで検討・発展させてきた、渦度―流れ関数法と圧力のポアソン方程式の連成による血管壁の拍動挙動解析手法について述べる。
２．

解析方法

2.1 解析モデル
本研究で用いた、直管を模した解析モデルを図 1 に示す。流体は均質な非圧縮性 Newton 流体、流れは軸対称流とし、流れ方向および半径方向の二次元で解析を行った。

図 1 解析モデル
2.2 基礎式
二次元円筒座標系における連続の式と Navier-Stokes 方程式をそれぞれ（１）及び（２）式に示す。
�����

���
��+ ��������+ ��+� �� �����
��(�� � ) = 0=−� ��� ��+��（１）� � ����+�� ���� ��+�� ���� ��

— 25 —CFD による血管壁拍動挙動における Vasomotion の再現解析�� ���+ ����� ���+ ���� �=−��� ��+��� ���� � ���+�� �� �� ��u+ : 軸方向速度m/s

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

ft+ : 時間sp:

圧力

Pa=N/m2ρ:

密度kg/m3ν:

動粘性係数m2/s+�� � �����− ����

（２）

なお、添え字の+は有次元量であることを表す。管内流れは（１）式を満たすことから、（３）及び（４）式を満たす Stokes の流れ関数�が存在する。

�� = −� ��

�� = −� ���=−

また、渦度�は（５）式で定義される。

（４）� ���� �
��

（３）

� ��

���

（５）

��

本研究では、流れ関数渦度法により流れ関数�と渦度�を計算することで流れ場の解析を行った。また、初期半径�� 、代表流速�� を用いて無次元化を行った。ここで、�� は（６）式に示す Reynolds 数の定義から、Re の値を定めることにより�� を決定した。

�� =

�� ��

（６）

�

以上により無次元変数及び無次元流れ関数�、無次元渦度�を導入することにより、（１）、（２）式を連立することで以下の Poisson 方程式及び渦度輸送方程式を以下の（７）及び（８）式のように導出できる。

����+���
��

ΩR =+���
��

�� �
��

−�

�

:

管出口:

管中心軸上:

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

��= 0,= 0,��������� ��− Ω=�

また、境界条件は以下のように定義される。
管入口� ��= 0,�� ��� ��+��� ��� ��� ��� ��（７）+��� ��� �+��� ���� �（８）=0= 0, � = �00����= 0,� = 0� � = 0,� = 0,

����

=0

� = 0� � = 0� � = −��0�定常流),� = −��0�� − ����)�拍動流),� =� �� � ��,=0� �� � ��

ここで、壁面上の流れ関数を余弦関数により振動させることで、拍動を再現している。

解析において、血流の拍動に伴う血管の拡張収縮により、血管壁は直線ではなく複雑な曲線形状を

— 26 —

SENAC Vol. 48, No. 4（2015. 10）持つ。そのような非直線形状を持った領域における流れを取り扱う場合、長方形格子を用いると曲線境界が階段形状に近似されるため、境界条件を正確に課すことは困難である。一方、境界に沿った格子を用いれば正確な境界条件を容易に課すことができる。そこで、座標変換により基礎方程式の変換を行い、解析領域を壁面形状で数値的に無次元化することにより物理領域である(R,Z)座標系から計算領域である(η,Z)座標系への変換による写像を施す。(R,Z)座標系における壁面形状を S(Z)で表すと、以下の（９）式で示される。

�=

�

�(�)

�� = �（９）これらを基礎方程式に適用すると、無次元基礎方程式である（７）及び（８）式は以下の（１０）及び（１１）式のように表される。

����

=

���

����� ��� ���� ���� ���� ����� ��� �

各係数は以下の通りである。

� ����
����������� =�� =�� �� ��� ����������� ���� ��� �= �� Ω� ��� Ω� − ���1��(�) ���1���
� �
�(�)�
��

��
��

（１０）− �������− ��� Ω （１１）

�

11
��(�)
� � �(�)
�−
�
2�
�
�
−
��(�)
�
�� �
�(�)�
�
��
�� = −

2� ��(�)�(�) ��

�� = 1

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

�� = −

2� ��(�)�(�) ���11��(�)� � �(�)
�� =
�
�
2��
�
−
��(�)
�
�� �
�(�)� ������ = − ���� =��� =1 �
�
��(�)

2 ��
��(�)� ��

��� = −2 ����(�)� ��2� ��(�) �� ���−����(�) �(�) �� �� ��（１０）、(１１)式を連立して�、�を計算することにより、Navier-Stokes 方程式の解を求める。なお、
�は SOR 法、�は陽解法により計算した。得られた�から、軸方向速度 u と半径方向速度 v は次のように求められる。

�=−���

��(�)� ��（１２）— 27 —CFD による血管壁拍動挙動における Vasomotion の再現解析

�=�

�−�

��(�)

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

��

+����

�

（１３）

次に、圧力の定式化について述べる。（２）式において圧力の半径及び軸方向の一次微分項を作り、（１４）式に示すポアソン方程式を求め圧力を計算した。
��������� ���� �� ������
��� = − � − �� � + 2+� � �− ������� ������∆は軸対称流における円筒座標系の微分演算子で、（１５）式で表される。

∆=���� �+� �� ��+���� �（１４）（１５）（１４）式における D+は、連続の式の値として以下のように定義する。

�� � � � � =�����+���+�� ���

（１６）

本来、流れ関数渦度法では連続の式が満たされるが、圧力式を差分化して解く際に離散化による誤差が発生する。より高精度な計算のため、D を含む項のうち計算に影響を及ぼすと考えられる時間微分項を式中に残してある[1]。（１４）式を無次元化及び写像、差分化した最終的な計算式を（１７）式に示す。

���� =0��� � ℎ�

1 + � � + �� ���(�)������� � ��(�) ������ + ������2� ��(�)� � � �(�) ������ − ����������+���−�� � ��ℎ������ �� �2ℎ����
2� ��(�) �������� − �������� − �������� + �������� ������ + �������
�−+4ℎ�ℎ�� �����
�−��+�11�������������
�����
�����
+ �+ �− ������2ℎ�ℎ�����− ������ � 2 ������ − ������1 �� ��(�) ������ − ������ ������ − ������ ��� � � ������� +�−+�����
2ℎ2ℎ2ℎ2ℎ�
� ���
�+���
�
���� �� ��(�) ������ − ������ ������ − ����������+ �−+� +� ���� ��2ℎ2ℎ������
（１７）

ここで、D の差分化において、t について離散化し、��れる圧力の値は、� ≠ 0であっても�������

= 0と置いている。このようにして求めら= 0となるように定めているため、境界条件や時間進行による誤差集積に対して D を常に小さな値に抑えておくことが出来る[1]。この手法は MAC 法で用いられるものであるが、本研究における解析手法は、流れ関数渦度法から得られる流れ関数から流速を計算し、時間発展させているという点で MAC 法とは異なっている。
先に述べたように、本解析手法では、流れ関数の壁面境界条件により流量及び流速に変動が生じ、それに伴い管内圧力が変動する。そこで、血管壁が管内圧力に応じて変形する効果を考える。血管壁の剛性はひずみに対して非線形性を示すことが知られており、血管壁の剛性を（１８）式に示す、fungによる理論式[9]より求めた。

�=�∗��∗ � ∗������� − �� �����(� − �∗ )�
�:血管壁に作用する応力[g/mm2]

:

ひずみ

（１８）

— 28 —

SENAC Vol. 48, No. 4（2015. 10）（１８）式をひずみで微分した（１９）式を計算することで、剛性 E を算出した。
���∗��= ∗ � �� � � � � �� � ��� �����(� � �∗ )�
（１９）
���
�
� � ∗�
�

ここで a、� ∗ 、�∗ は定数であり、体重 15~25kg のイヌの胸部大動脈を用いた引張強度試験による実験結果[10]から、以下の値を用いた。

a�∗

�∗

:

1.9:25.83g/mm2:

1.6

内圧 p、壁厚 h、半径 R の薄肉円管モデルにおける円周応力は以下のようになる。

�=���

（２０）

また、円周方向のひずみ�と応力�の関係は以下のようになる。

�� = ���（２１）以上より、半径は（２２）式で表される[11]。なお、以下のように漸化式の形とすることで、圧力変動に対し各ステップの半径を計算することができる。

��ℎ���� =�� = ���� � ���������� �:��:

管壁の厚さ m:

圧力の増分 Pa=N/m2（２２）

管壁の剛性増分 Pa=N/m2

壁面形状関数として定義された S(Z)は（２２）式により計算される。これにより、各軸方向位置における半径を用いて写像を行うことで、計算領域を長方形格子に変換している。
次に、解析のフローチャートを図 2 に示す。第一節で述べたように、本解析では渦度流れ関数法と圧力のポアソン方程式を連成して解く手法を取っている。一般的に、レギュラー格子を用いた差分法
による流体解析では、求めた圧力から流速を時間発展させることにより連続の式が完全には満足されないという問題がある。しかし、本解析手法では流れ関数渦度法により時間発展させるため、レギュラー格子を用いながらも連続の式を十分に満たし、かつ圧力を計算することが可能となる。

— 29 —CFD による血管壁拍動挙動における Vasomotion の再現解析

図 2 フローチャート
2.3 解析条件
2.3.1

Vasomotion のモデル化

血管運動には、心拍動の周期とは異なる、自律神経活動に起因する自然発生的で律動的な長周期的収縮拡張運動(Vasomotion)が存在する[12]。一般に、Vasomotion の生理学的意義については多くの研究がなされている[13]。しかし、Vasomotion が血管壁拍動挙動に及ぼす影響についてはまだ解明されていない。そこで、Vasomotion が関与する血管壁運動を数値解析において再現するため、拍動の周波数�� より低い周波数�� で振動する関数を用いて、（２３）式で示す拍動条件を考えた。
��Ψ���� = −1.0 × (� − ����� �) × �

��−�������� ���

��� = �� ��（２３）

この拍動条件を用いて、Vasomotion が関与する血管壁の拍動挙動に与える影響について検証した。
2.3.2 粘弾性モデル
血管壁は弾性的性質が劣化すると、その因子であるコラーゲン等が有する粘弾性的性質により、クリープや応力緩和に起因する粘弾性的特性が発現する。この時、血管壁がどのような拍動挙動を示すかを明らかにすることは血管疾患の診断において重要である。本解析を行うために、血管壁の応力ひ— 30 —SENAC Vol. 48, No. 4（2015. 10）

ずみ構成則として図 3 に示すような粘弾性ばねモデルを構築した。このモデルは、粘弾性 Voigt モデ
ルと Maxwell モデルを組み合わせたものであり、モデル血管における血管壁のクリープ特性と応力緩和特性を同時に再現している[5]。この粘弾性モデルに引張応力を負荷した場合の構成式は、（２４）、（２５）式で表される。
�σ = �� � � �� �� � ��� �� � ���E=�������= �� � �� ��� �� � ���� � ����

��

:

弾性係数:

粘性係数:

ひずみ速度（２４）（２５）

図 3 粘弾性ばねモデル
ところで、粘弾性特性には、応力負荷時と除荷時で応力ひずみ曲線が異なる経路をたどる、ヒステリシス特性が存在する。より詳細に血管壁粘弾性をシミュレートするため、ヒステリシス特性を再現する条件を導入した。除荷条件として、最大ひずみ���� からの減少量���� を導入し、（２４）式をΔσとして最大負荷応力���� からの減少量と定めた。これにより、負荷曲線を 180°回転した形の除荷曲線が決定し、ヒステリシス特性が再現される。
2.4 カオス理論に基づく不安定挙動解析

過去の研究において、血管壁粘弾性の発現により血管挙動の不安定性が増加することが臨床において示されている[4]。著者の一人が提案してきた DPC 法[2-5]では、測定した血管壁ひずみ速度から動脈硬化進行度を定量的に評価している。さらに、血管壁挙動の振動成分をより詳細に解析するために、アトラクタ解析による軌道不安定性の定量化が行われている[4]。図 4～6 に示すように、DPC 出力波形に対しアトラクタ解析を行うことで、血管壁挙動の拍動軌跡が得られる。また、以下の式で示されるエントロピーS を計算することで、軌道の不安定性を定量化することが出来る。

— 31 —CFD による血管壁拍動挙動における Vasomotion の再現解析(a) 健常例(b)動脈硬化症例

図 4 超音波ドップラー法で検出された血管壁拡張時の血管壁変形速度非侵襲測定例[2]

(a)

(b)

V(t+⊿t)V(t+⊿t)

図 5 拍動下の(a)血管壁速度非侵襲測定結果と(b)拍動軌跡特性[4]

V(t)

V(t)

(a) 健常例(b)動脈硬化症例
図 6 血管壁の拍動軌跡[4]
�

�� :

� :

� � � � �� �����
���

i 番目の領域に軌道が落ちる確率
領域分割数

（２６）

— 32 —

SENAC Vol. 48, No. 4（2015. 10）

健常者において、血管壁の拍動挙動には図 6(a)に示されるように、サーキットループ(リミットサイクル)の軌跡に長周期的変動が重畳されていることが分かる。動脈硬化症の血管においては、このループが半楕円状になり、また拍動軌道に乱れが生じていることが分かる。また、血管壁拍動軌道の乱れがエントロピーの増加として発現している [4] 。本稿では、数値解析結果に対し血管壁拍動軌跡を解析し、臨床結果との比較を行った。
３．

解析結果

図 7～9 に弾性変形モデルと Vasomotion を考慮した弾性変形モデル、粘弾性モデルのそれぞれについて数値解析を行い、得られた血管壁ひずみ速度波形及び拍動軌跡を示す。また、これらの各解析条件における、血管壁ひずみ速度波形から求めたエントロピーS を表１に示す。
弾性変形モデルでは、血管壁ひずみ速度波形は健常例のような正弦波の形を示し、また波形の経時的変化も見られなかった。このような定拍動振幅条件においては、拍動軌道は図 7(b)のようにただ一つの軌道を示すことが解析により示された。
Vasomotion を考慮した弾性変形モデルでは、通常の弾性変形モデルと比べ血管壁ひずみ速度波形の形状的変化は見られなかったが、時間進行によりその振幅が異なる波形が得られた。この速度波形に対する拍動軌跡には、図 8(b)に示されるような長周期的なゆらぎが確認された。このように、Vasomotion の導入により拍動軌跡に長周期的軌道が重畳され、図 6(a)で示すような健常者の臨床例に一致することが分かった。すなわち、健常血管においても血管壁拍動挙動に Vasomotion が存在することが CFD により検証された。
次に、粘弾性効果を導入したモデルでは、図 9 に示すように時間進行に伴う血管壁ひずみ速度波形の拍動振幅に変動がみられ、また拍動軌道のゆらぎが確認された。このモデルのように、定拍動振幅条件においても粘弾性の効果により粘弾性特有の Vasomotion が発現することが解析により示された。
さらに、Vasomotion を考慮した粘弾性モデルでの解析結果を図 10 に示す。この解析条件では、通常の Vasomotion によるゆらぎに対し粘弾性に起因する Vasomotion が重畳され、より顕著な軌道の乱れが確認された。また表 1 から、粘弾性の発現による軌道の乱れがエントロピーの増加として反映されていることが示された。
４．

結言

渦度－流れ関数法と圧力のポアソン方程式を連成し血管壁の拍動挙動を解析する手法により、健常な血管壁拍動における長周期的振動の存在と、壁粘弾性発現による特有の長周期的振動という二つのVasomotion の存在を再現できた。今後、さらに本手法を拡張発展させることで、詳細な血管疾患との対応が可能となるだろう。

(10-5)V(t+⊿t)Wallspeed

3210500000

Timestep

600000(a) 血管壁ひずみ速度波形
図 7 弾性変形モデル

V(t)(b) 拍動軌跡

— 33 —CFD による血管壁拍動挙動における Vasomotion の再現解析(10-5)4V(t+⊿t)

Wallspeed

3210500000

600000TimestepV(t)(a) 血管壁ひずみ速度波形
図8(b) 拍動軌跡

Vasomotion を考慮した弾性変形モデル(10-5)2V(t+⊿t)

Wallspeed

310500000

600000TimestepV(t)(a) 血管壁ひずみ速度波形(b) 拍動軌跡

図 9 粘弾性モデル(10-5)V(t+⊿t)Wallspeed

3210500000

Timestep

600000V(t)(a) 血管壁ひずみ速度波形
図 10(b) 拍動軌跡

Vasomotion を考慮した粘弾性モデル
表 1 エントロピー
解析モデル

エントロピー

弾性変形モデル

4.417

Vasomotion を考慮した弾性変形モデル

4.942

粘弾性モデル

4.817

Vasomotion を考慮した粘弾性モデル

5.186

謝辞
本解析手法の開発にあたり、東北大学サイバーサイエンスセンターの SX シリーズを使用させて頂いたことを感謝する。

— 34 —

SENAC Vol. 48, No. 4（2015. 10）

参考文献[1] 河村哲也

流体解析１、1996,朝倉書店[2] A.T.Yokobori,Jr., T.Ohkuma, S.Sasak, H.Yoshinari, T.Yokobori, H.Ohuchi and S.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. 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 andH.Yamashita, Bio-Medical Materials and Engineering, 23, (2013), pp.75−91.
[6] A.T.Yokobori.Jr., Y.Nishikawa, The analysis of pulsatile behavior of blood vessel wall withaneurysm based on hybrid analysis of computer fluid dynamics (CFD) and solid mechanics,
Proc. of the Ann Meeting of Atherosclerosis and Thrombosis, (2013)[7] A.T.Yokobori.Jr., Y.Nishikawa, The analyses of mechanical behaviors of blood vessel wall withaneurysm at the progressive stage based on simulate experiment and computer fluiddynamics(CFD), Proc. of the Ann Meeting of Atherosclerosis and Thrombosis, (2014)[8] A.T.Yokobori.Jr., K.Kudo, Long Periodic Fluctuation of Blood Vessel Wall Caused by theViscoelasticity Under Pulsatile pressure Flow, Proc. of the Ann Meeting of Atherosclerosis
and Thrombosis, (2015)[9] Y.C.Fung, Appl. Mech. Rev., 1968,21(1),pp1-20.
[10] 横堀壽光，大熊恒郎，前山俊秀，横堀武夫，大内博，奈良茂樹，佐々木久雄，葛西森夫，生体材料，1986,4,pp27-32.
[11] 岡小天，ながれ，1993，12，pp249-254.
[12] 南山求，中野厚史，日本レーザー医学会誌，1994,9,pp27-38.
[13] T.Sakurai, N,Terui, Am J Physiol Heart Circ Physiol 291, (2006), pp.1761-1767