第 25 回数値流体力学シンポジウムD04-2

赤血球を多数含んだ微小血管での血流解析
Numerical analysis of blood flow including multiple red blood cells in capillary vessel
伊井 仁志, 東大, 文京区本郷 7-3-1, sii@fel.t.u-tokyo.ac.jp:
杉山 和靖, 東大, 文京区本郷 7-3-1, sugiyama@fel.t.u-tokyo.ac.jp:
高木 周, 東大/理研, 文京区本郷 7-3-1/和光市広沢 2-1, takagi@mech.t.u-tokyo.ac.jp:
松本 洋一郎, 東大, 文京区本郷 7-3-1, ymats@fel.t.u-tokyo.ac.jp:

Satoshi Ii, The University of Tokyo, 7-3-1 Hongo Bunkyo-ku
Kazuyasu Sugiyama, The University of Tokyo, 7-3-1 Hongo Bunkyo-kuShu Takagi, The University of Tokyo/RIKEN, 7-3-1 Hongo Bunkyo-ku/2-1 Hirosawa Wako
Yoichiro Matsumoto, The University of Tokyo, 7-3-1 Hongo Bunkyo-ku
A numerical simulation is conducted for a blood ﬂow in a capillary vessel, including not only multiple red bloodcells (RBCs) but also platelets. The recently-developed full Eulerian approaches for coupling the ﬂuid with both
elastic structure and membrane are employed for dealing with the RBC, platelet and vessel wall. A pressureinduced periodic ﬂow is imposed, and then the numerical run is carried out until a certain time, where the ﬂow
is developed well. The numerical results indicate that the platelet motions are strongly aﬀected by a dynamic
motion of a RBC and its multiple behaviors as well.

はじめにここ数十年，血流 (あるいは血漿) と赤血球・血小板・白血球・血管壁などの弾性物との力学的な連成を取り扱う様々な解析手法が提案されている．中でも，Peskin(17)によって提案された Immersed Boundary(IB) 法では，流体中を大変形かつ大移動する弾性物の取扱いが容易であり，その適用性の高さから現在でも幅広く用いられている．IB 法の発想では，移動・変形する弾性物は物質移動点で表現され，そこで得られる特異力をオイラー記述された固定格子に分配することで流体と弾性物の連成が行われる．物質点を陽に保持しているため，数値拡散の抑制や初期配置に対する相対的な変形量が容易に分かるという利点がある一方，弾性物が極端に変形した場合，固定格子と物質移動点との空間解像度比による数値的な不安定性や，物質移動点から構成されるラグランジュ格子の破綻など，これらを防ぐために格子再構築などの安定化が必要となり，それに伴う誤差の蓄積が問題となる場合がある．
このような問題を避けるため，固定格子のみを用いた完全オイラー型定式化による連成解析手法が提案されている (15,5,2,3) ．完全オイラー型手法では，物体を物質点の集合としてラグランジュ的に表現するかわりに，固定格子において導入された媒質量によりオイラー的に表現する．
この際，変形量を記述する量も導入し，これらの時間発展をオイラー的に求める．近年，生体をターゲットとした完全オイラー型定式化による流体と超弾生体 (22,16,23,11,20) ，また流体と超弾性膜に対する連成解析手法 (12,20) が提案されており，界面を陽に記述する従来手法に比べ本質的に数値拡散が潜在するものの，格子解像度を上げるにつれ従来手法に収束していることが確認されている．完全
オイラー型手法では，前述した格子再構築に関する問題が存在しない他に，全ての変数を固定格子において定義できるという利点から，大規模計算における並列化効率の向上にも期待されている．
そこで本研究では，完全オイラー型手法を多数の赤血球・血小板を含む微小血管での血流解析に適用し，その適用性を調べるとともに，赤血球と血小板の力学的な挙動に関して考察を行う．

1.

2.

解析手法

2.1 オイラー型定式化による連成モデル
混合近似された，いわゆる 1 流体モデルの定式化において，流体，弾性体また弾性膜内部の各媒質相をそれぞれ上添字 (f ) : f luid，(s) : solid および (m) : membrane1と表すと，連成モデルは以下のように記述される．

∂φ(α)+ (v · ∇)φ(α) = 0, (α = s, m),∂t()φ(ρ(f )=1−φ∂v+ (v · ∇)v∂t)

(s)

−φ(m)

(1),

∇ · v = 0,= −∇p + ∇ · τ + Fs −(2)∆Pex , (3)L∂B(s)+ (v · ∇)B(s) = B(s) · ∇v + ∇vT · B(s) ,∂t(4)

(m)

∂Bs
∂t+ (v · ∇)B(m)= B(m)· ∇s v + ∇s vT · B(m)sss , (5)(m)

∂Js∂t+ (v · ∇)Js(m) = Js(m) ∇s · v

(6)

(m)

∂κs0(m)+ (v · ∇)κs0 = 0.
∂t(7)ここで，φ(α) は媒質 α(= f, s, m) の体積占有率 (VOF関数)(9) ，v は速度ベクトル，p は不定圧力，τ は応力テン
ソル，Fs は表面特異力ベクトル，−∆P/L は駆動圧力勾配，B(s) は媒質 s に関係するバルクでの左 Cauchy-Green(m)(m)(m)変形テンソル，また Bs ，Js および κs0 はそれぞれ媒質 m に関係する表面左 Cauchy-Green 変形テンソル，初期配置と現配置との表面積比，また初期配置での平均曲率である．式 (5) および (6) において，∇s は射影面での微分演算子であり，媒質界面の単位法線ベクトル n から構成される表面射影テンソル P = I − nn を用いると，∇s = P · ∇ となる．ここで，媒質界面を陽に保持していないため，単位法線ベクトルは n = ∇φ/|∇φ| より求める．各変数の詳細な定義は，例えば，文献 (21,1,19) に，またより詳細な定式化は文献 (23,12,20) に示されている．
先行研究 (12) に対し，本研究では媒質の表面積比Js を√導入している．本来，Js = (tr(Bs )2 − tr(B2s ))/2 の定義よりその時間発展を考慮する必要はないが，定義における積の項が移流計算などの数値的な影響により不安定

Copyright c 2011 by JSFM

第 25 回数値流体力学シンポジウム
D04-2

解析結果
Fig.1 に示すよう，直方体の計算領域に円筒管を考え，その中に赤血球を模擬した弾性膜，また血小板を模擬した弾性体をランダムに配置する．ここで直方体と円筒管の間に血管壁を模擬した弾性体を考慮している．圧力勾配を体積力として主流方向に作用させ系を駆動する．主流方向 x には周期境界条件，また y および z の境界条件には固定壁を課す．用いた計算条件および物性値は Table
1 にまとめる．この条件において，赤血球の容積率 (ヘマトクリット) は約 20%となっている．本解析では血小板・血管壁と 2 種類の弾性体を定義しているため，各弾性体に対し φ および B を導入する．また血小板と赤血球に関して，数値的な合体を防ぐため各個体に対して φ を導入する．ただし，B や Bs は各媒質で 1 種類とする．

性を引き起こしている可能性がいくつかの数値テストから確認されているため，Js を導入しその時間発展を式 (6)から求めている．
1 流体モデルでは，密度および応力テンソルは以下のように VOF 関数で混合化される．

3.

ρ = φ(f ) ρ(f ) + φ(s) ρ(s) + φ(m) ρ(m) ,(8)τ = φ(f ) τ (f ) + φ(s) τ (s) + φ(m) τ (m) .

また，VOF 関数が計算格子幅で滑らかに遷移していると仮定し，界面法線方向のデルタ関数を δn ≈ |∇φ| と近似
すると，表面特異力ベクトル Fs は媒質 m の膜応力ベク(m)トル fs を用い，以下のように与えられる．

Fs = |∇φ|(m) fs(m)(9)

2.2 構成方程式
流体はニュートン流体，また弾性体は Neo-hookean モデルおよびニュートン粘性則に従うと仮定すると，応力テンソルは以下のように与えられる．
τ (f ) = 2µ(f ) D,
τ (s) = 2µ(s) D + G(s) (B(s) − I),(10)τ (m) = 2µ(m) D,(α)

Fig. 1: An initial conﬁguration of the RBCs, platelets
and vessel wall.

(s)

ここで µ
(α = f, s, m) は各媒質の粘性係数，G は弾性体の横弾性係数，また D = (∇v + ∇vT )/2 はひず
み速度テンソルである．
膜応力は赤血球の変形挙動を表すのに良く用いられるEvans-Skalak モデル (6) と曲げ剛性モデル (18) を導入する．ここで，上添字 (m) は簡単化のため省略する．

fs = ∇s · τs + qn,Fig.2 に結果を示す．赤血球・血小板が変形しながら流れにのって移動しており，この時，赤血球は実際の血流で観察されている，いわゆるパラシュート形状やスリッパ形状 (7) の様相をとっていることが確認できる．時間進行するにつれ，流体力学的な効果により赤血球集団が軸方向に集まっており，血管壁との間に cell-free layer と呼ばれる血漿相が形成されている．それに伴いサイズの小さい血小板の挙動が影響を受けている．そこで，血小板重心の血管動径方向の時間履歴を調べたところ (Fig.3)，
動径位置 4 ≤ r ≤ 7 付近では血小板の存在する機会が非常に少なくなっている．各時間における赤血球の動径方向の空間占有率を Fig.4 に示す．流れが発展するにつれ
cell-free layer が大きくなり，軸中心付近の赤血球の占有率が増加している．これらの事実より，血管壁付近の血小板は赤血球集団の高密度な空間占有によるバリアを越えられず cell-free layer に捕獲される．一方，軸中心付近の血小板は動径方向の移動が赤血球の高密度相により制限される．これらは二次元における解析結果からも支持される (4) ．

(11)と膜応力が与えられ，()
Λ2 + 1
CsG+E
Λ
−CP,
ss 1s
(Λ1 + 1)2
Λ1 + 1tr(Bs )− 1,
Λ1 = Js − 1, Λ2 =2Js
τs =(12)(Gs = P · Bs · P),
また，

q = (∇s · m) · P,

まとめ
本研究では，オイラー型連成解析手法を微小血管中の血流解析に適用し，赤血球や血小板の流動挙動を解析した．結果より，赤血球の集団挙動が血小板挙動に強く影響を与えることが分かった．その際，血管軸付近に赤血球の高密度な空間占有相が形成され，それが動径移動のバリアとして働くことにより，血小板の存在位置あるいは領域が 2 分化されていることが示唆された．今後はサンプル数を増やして統計的な向上を行い，より詳細な検証を行っていく．

4.

(13)

m = Eb (∇s n + κs0 P) ,である．Es ，Cs はそれぞれ膜弾性係数と面積膨張係数，また Eb は曲げ剛性係数である．
なお式 (11) において，τs = P · τs の関係より，実装では ∇ · τs として評価している．

2.3 数値モデル
速度と圧力のカップリングにはプロジェクション法 (8)を用いる．数値的に ∇ · v = 0 を満たす速度場により，式(1) および式 (4)-(7) の移流項をそれぞれ MTHINC 法 (13)
また五次精度 WENO 法 (14) ，また右辺項を二次精度中心差分により離散化し，二次精度ルンゲ・クッタ法より時間発展を行う．その後，得られた物理量を用い，運動方程式 (3) を二次精度中心差分法で離散化する．ここで，移流項にはオイラー陽解法，また粘性項にはオイラー陰解法を適用する．
2Copyright c 2011 by JSFM

第 25 回数値流体力学シンポジウム
D04-2

0.6
t=0.0 [ms]t=7.5 [ms]t=30.0 [ms]t=75.0 [ms]volume fraction [-]

0.5

0.4

0.3

0.2

0.1

0012

3456

78910

radial position [µm]Fig. 4: Volume fraction versus axial position of themultiple RBCs at diﬀerent time.

Tab. 1: Computational conditions.

Fig. 2: The snapshots of the numerical results att = 7.5, 39 and 67.5 [ms].

109

radial position [µm]8765432
10010203040506070t [ms]Fig. 3: Developments of the axial positions for the respective platelets.

3Domain [µm3 ][0,44]×[0,22]×[0,22]

Pressure gradient [kPa/m]
Spatial resolution [-]-400
320×160×160

Time increment [µs](Plasma)density [kg/m3 ]1.5

viscosity [cP](Wall)

1.2

radius [µm]density [kg/m3 ]viscosity [cP]

1010001.2

elastic modulus [N/m2 ](Platelet)

5num. of pieces [-]radius [µm]density [kg/m3 ]1011000

viscosity [cP]elastic modulus [N/m2 ]
(RBC)

1.25

number of pieces [-]radius [µm]30
3.9

inner liquiddensity [kg/m3 ]viscosity [cP]1000
6membrane (10,19)
elastic modulus [dyn/cm]

0.006

dialation modulus [dyn/cm]
bending modulus [dyn·cm]

0.6
2×10−121000

Copyright c 2011 by JSFM

第 25 回数値流体力学シンポジウム
D04-2

参考文献
(1) D. Barthés-Biesel and J.M. Rallison, The timedependent deformation of a capsule freely suspended in a linear shear ﬂow, J. Fluid. Mech., 113(1981) 251-267.

(17) C.S. Peskin, Flow patterns around heart valves: anumerical method, J. Comput. Phys., 10 (1972)
252-271.
(18) C. Pozrikidis, Eﬀect of bending stiﬀness on the deformation of liquid capsules in simple shear ﬂow, J.
Fluid. Mech., 440 (2001) 269-291.

(2) G.H. Cottet and E. Maitre, A level set method for
ﬂuid-structure interactions with immersed surfaces,
Math. Model. Meth. Appl. Sci., 16 (2006) 415438.

(19) C. Pozrikidis, Modeling and Simulations of Capsules and Biological Cells, Chapman & Hall, CRC:Boca Raton, 2003.

(3) G.H. Cottet, E. Maitre, T. Milcent, Eulerian formulation and level set models for incompressible
ﬂuid-structure interaction, Math. Modeling andNumer. Anal., 42 (2008) 471-492.

(20) S. Takagi, K. Sugiyama, S. Ii, Y. Matsumoto, A Review of Full Eulerian Methods for Fluid StructureInteraction Problems, J. Appl. Mech., accepted(2011).

(4) L.M. Crowl and A.L. Fogelson, Computationalmodel of whole blood exhibiting lateral platelet motion induced by red blood cells, Int. J. Numer.
Meth. Biomed. Eng., 26 (2010) 471-487.

(21) R. Skalak, A. Tözeren, P.R. Zarda and S. Chien,
Strain energy funcion of red blood cell membranes,Biophys. J., 13 (1973) 245-264.

(5) T. Dunne, An Eulerian approach to ﬂuid-structureinteraction and goal-oriented mesh adaptation, Int.
J. Numer. Meth. Fluids, 51 (2006) 1017-1039.

(22) K. Sugiyama, S. Ii, S. Takeuchi, S. Takagi and Y.
Matsumoto, Full Eulerian simulations of biconcave
neo-Hookean particles in a Poiseuille ﬂow, Comput.
Mech., 46 (2010) 147-157.

(6) E.A. Evans and R. Skalak, Mechanics and Thermodynamics of Biomembranes, CRC, Boca Raton,
FL, 1980.

(23) K. Sugiyama, S. Ii, S. Takeuchi, S. Takagi, Y. Matsumoto, A full Eulerian ﬁnite diﬀerence approachfor solving ﬂuid-structure coupling problems, J.
Comput. Phys., 230 (2011) 596-627.

(7) P. Gaehtgens, C. Duhrssen and K.H. Albrecht, Motion, deformation, and interaction of blood cells
and plasma during ﬂow through narrow capillarytubes, Blood Cells, 6 (1980) 799-817.
(8) F.H. Harlow and J.E. Welch, Numerical calculation of time-dependent viscous incompressible ﬂow
of ﬂuid with free surface, Phys. Fluids, 8 (1965)
1834-1845.
(9) C.W. Hirt, B.D. Nichols, Volume of ﬂuid (VOF)
methods for the dynamics of free boundaries, J.
Comput. Phys., 39 (1981) 201-225.
(10) R. Hochmuth, R. Waugh, Erythrocyte membraneelasticity and viscosity, Annual Review of Physiology, 49 (1987) 209-219.
(11) S. Ii, K. Sugiyama, S. Takeuchi, S. Takagi and Y.
Matsumoto, An implicit full Eulerian method forthe ﬂuid-structure interaction problem, Int. J. Numer. Meth. Fluids, 65 (2011) 150-165.
(12) S. Ii, X. Gong, K. Sugiyama, J. Wu, H. Huang,S. Takagi, A full Eulerian ﬂuid-membrane coupling
method with a smoothed volume-of-ﬂuid approach,Commun. Comput. Phys., accepted (2011).
(13) S. Ii, K. Sugiyama, S. Takeuchi, S. Takagi, Y. Matsumoto, F.Xiao, An interface capturing methodwith a continuous function: the THINC methodwith multi-dimensional reconstruction, J. Comput.
Phys., revised (2011).
(14) G.S. Jiang, C.W. Shu, Eﬃcient implementation ofWENO schemes, J. Comput. Phys., 126 (1996)
202-228.
(15) C. Liu, N.J. Walkington, An Eulerian description ofﬂuids containing visco-elastic particles, Arch. Rational Mech. Anal., 159 (2001) 229-252.
(16) N. Nagano, K. Sugiyama, S. Takeuchi, S. Ii, S.
Takagi and Y. Matsumoto, Full-Eulerian FiniteDiﬀerence Simulation of Fluid Flow in Hyperelastic
Wavy Channel, Journal of Fluid Science and Technology, 5 (2010) 475-490.
4Copyright c 2011 by JSFM