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

完全オイラー型定式化による流体・弾性膜連成問題の解析
A full Eulerian approach for the fluid-elastic membrane coupling problem
伊井 仁志, 東大工, 文京区本郷 7-3-1, sii@fel.t.u-tokyo.ac.jp:
杉山 和靖, 東大工, 文京区本郷 7-3-1, sugiyama@fel.t.u-tokyo.ac.jp:

GONG Xiaobo, SJTU, 800 Dongchuan Road, Minhang District, Shanghai 200240, China, x.gong@sjtu.edu.cn:
高木 周, 東大工/理研, 文京区本郷 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-ku
Xiaobo GONG, Shanghai Jiaotong University, 800 Dongchuan Road, Minhang District, Shanghai 200240, ChinaShu 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 full Eulerian fluid-elastic membrane coupling solver is proposed. In order to describe the multi-component
materials, the volume-of-fluid approach is introduced. Furthermore, the surface deformation gradient tensor isupdated to deal with the membrane deformation in the Eulerian fixed-grid. The present method is validated fora problem in which a membrane capsule is deformed in shear flow, and is applied for a pressure-driven flow withthe biconcave membrane capsules.

1.

はじめに

2.

計算モデル

2.1 弾性膜モデル (11),(1),(10)膜弾性による表面張力のコーシー応力 τs は，表面ひずみエネルギー Ws が表面の変形量で決定されることより，以下のように記述される．

生体内において，血流における赤血球・血小板などの弾性膜と流体との連成問題の解析は重要な課題である．通常，そのような流体と弾性膜の連成解析手法では Peskin
によって提案されている Immersed Boundary 法 (8) が広く用いられる．Immersed Boundary 法では，物質点の集合として陽に記述された界面の表面張力を，オイラー場で解かれる流体の運動方程式に作用させ，連成させている．物質点を陽に保持しているため，初期配置からの変形の相対関係より，膜面応力を直接的に表現できるが，界面が極端に変形した場合，流体場での格子点に比べ界面の物質点数の密度が変わってしまい，格子解像度以下の高波数の変形モードの数値不安定性をさけるため，メッシュ再構築などの安定化が必要となる．

√ ∂Ws2 ∂Ws
Bs + 2 c2P.
τs = √c2 ∂c1∂c2(1)

ここで，P は界面法線ベクトル n で表される表面投影テンソル P = I − nn，また，Bs は表面で定義された左コー
シー・グリーン変形テンソルであり，c1 , c2 はその不変量を用いて，c1 = tr(Bs ) − 1, c2 =一方，流体とバルク超弾性体の完全オイラー定式化による連成解析手法が近年提案されている (12,13,5,7) ．完全オイ
ラー型手法では，物体を物質点の集合としてラグランジュ的に表現するかわりに，多媒質をオイラー格子で表現するために媒質の格子占有率 (Volume-Of-Fluid/VOlumeFraction:VOF 関数)(4) を用い，また変形量を記述する左
コーシー・グリーン変形テンソルを導入し，超弾生体と流体との連成解析に成功している．界面を陽に記述する
従来手法と比べ，若干，拡散的な振る舞いを示すが，格子解像度を上げるにつれ従来手法に収束することが確認されている．完全オイラー型手法では，本質的に格子生成および物質的の再構築のテクニックが不要であり，また，全ての変数がオイラー固定格子で時間発展が求められるため，将来的に大規模計算を行った際の並列化効率の観点などから，非常に有用な手段として考えられる．

)
1(tr(Bs )2 − tr(B2s ) ,2(2)と定義している．表面ひずみエネルギーのモデルには，
Neo-Hookean モデル (1) ，()Es1
Ws =c1 +,(3)6c2と赤血球で用いられる Skalak モデル (11) ，Ws =)Es ( 2c1 + αc22 − 2(α + 1)c2 + α − 1 ,8(4)などが提案されている．ここで，Es は膜弾性係数，α はSkalak モデルにおける膜膨張に関する係数であり，通常
α  1 とされる．以下のように定義される表面変形勾配
Fs より，そこで，本研究では完全オイラー型手法による流体と弾性膜の連成解析手法を提案する．膜の内外を区別するために VOF 関数を用い，膜変形をオイラー場で扱うために変形テンソルを導入する．また，VOF 関数の移流にMulti-dimensional THINC 法 (6) を用いる事で，界面の時間変化に伴う数値拡散を抑えつつ，滑らかに変化する
VOF 関数の特性を利用する．検証として，せん断場中におかれた球形カプセルの変形挙動，また圧力勾配で駆動される流れにおける赤血球を模擬した弾性膜の時間挙動を解析し，考察を行う．

Fs = P · F · PR = P ·∂x· PR ,∂X(5)

表面左コーシー・グリーン変形テンソル Bs は，

Bs = Fs · FTs = P · F · PR · FT · P,

(6)となり，変形勾配 F と初期配置と現在配置の表面投影テ
ンソル PR , P を用いて表される．

1

Copyright c 2010 by JSFM

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

また，曲げ剛性に関するモデルも提案されており (10) ，そのコーシー応力テンソルを qn とすると，q は曲げ剛性係数 Eb と初期配置での平均曲率 κRm を用い以下のようになる．

q = (∇s · m) · P,()
m = Eb ∇s n + κRmP .

Bs = P · G · P,

(7)

∂G+ (v · ∇)G = ∇vT · G + G · ∇v,∂t(9)∇ · v = 0,)(∂v+ (v · ∇)vρ∂t()
= −∇p + ∇ · µ(∇v + ∇vT ) + f .

(10)

3.

(11)
Re =

本研究では，膜にかかる特異な表面張力は，近似デルタ関数として δ ≈ |∇φ| を用い，

ργ̇d2
µγ̇d2U, Ca =, γ̇ =.
µEsHz

(19)

ここで，U は x 方向の速度，また Hz は鉛直方向の高さである．本計算では，ρ = 1，Re = 0.01，γ̇ = 1，また，
Ca = 0.0125 × 2n とし，n=0,1,2,3,4 を用い解析を行う．
計算格子 Nx = 128，Ny = 64，Nz = 128 における定常状態の計算結果を Fig.1 に示す．粘性はレイノルズ数を固定しているので，キャピラリー数を大きくするにつれ膜が柔らかくなり，その挙動が捉えられている．ここで，y 方向の中心断面における形状を楕円適合させ，長径l と短径 a で表される Taylor の変形パラメータ (14) D =(l−a)/(l+a) の時間履歴を Fig.2 に示す．また，Immersed
Boundary 法 (2) の結果，及び，Stokes 流体での Boundary
Element 法 (9) の結果を示す．これらの手法では，界面を陽に記述しているため，提案した完全オイラー型定式化より精度が良いと予想される．オイラー型手法では，界面応力が拡散的に表現されているため，膜が硬くなるにつれ，若干，本来より柔らかく表現される傾向にあるが，格子解像度を上げるにつれ，従来手法の解に収束する様子が確認できる．

(12)と評価する．ここで，曲げ応力に関して，初期配置での平均曲率 κRm を含むので，(13)Rと変形した評価式を用いる．κ̃Rm = |∇φ|κm であり，数値R拡散を抑えるために，κ̃m を予報変数として求める．よって，時間発展式は，

κ̃R
mR+ (v · ∇)κ̃Rm = −κ̃m n · ∇v · n,∂t

計算結果

3.1 せん断場中における弾性膜計算領域 Hx × Hy × Hz = [−2, 2] × [−1, 1] × [−2, 2]
において，鉛直方向の上下境界に速度を与えたせん断場中に，直径 d = 1 の球形の弾性膜 (Neo-Hookean モデル)を配置し，その挙動を検証する．ここで，無次元数としてレイノルズ数 Re，キャピラリー数 Ca，せん断率 γ̇ を以下のように定義する．

ここで，φ は媒質の体積平均値 (VOF 関数)，v は速度ベクトル，また，p は不定圧力である．混合モデルでは，各相の密度 ρ1 , ρ2 ，粘性係数 µ1 , µ2 は VOF 関数を用いて以下のように平均化される．

q̃ = (∇ · m̃) · P,()m̃ = Eb |∇φ|∇n + κ̃RmP .

(18)

2.3 計算手法
本計算では，速度と圧力のカップリングにプロジェクション法 (3) を用いる．VOF 関数の移流計算手法に界面捕獲手法の一種である Multi-dimensional THINC 法 (6) ，また G の移流項には三次精度風上型有限体積法，他の項には二次精度中心差分を用いる．時間離散には，移流項は二次精度 Adams-Bashforth 法，粘性項には CrankNicolson 法を用いる．流体では変形に対する拘束がない
ため G は数値的に発散する可能性があるので，|∇φ| ≤ の判定により，G = P と再初期化する． は問題に依存し，本研究では  ≈ 0.001 − 0.1 ほどとしている．

ここで，Γ は界面，δΓ は三次元デルタ関数，また，H は各相を表す特性関数である．
本研究では，混合近似された基礎方程式を用いる．

∂φ+ (v · ∇)φ = 0,∂t(17)となる．

Γf = |∇φ|∇s · (τs + qn) = ∇ · (|∇φ|τs + q̃n),(16)と表面の左コーシー・グリーン変形テンソルを計算する．
よって，G の時間発展式は，Ḟ = ∇vT · F の関係式を用い，2.2 オイラー型定式化運動方程式において，表面張力はデルタ関数を用いて体積力 f として作用される．
IδΓ fs dΓ = |∇H|fs .
(8)
f=µ = µ1 φ + µ2 (1 − φ).

G ≡ F · PR · FT ,と定義し，G の時間発展を求め，

ここで，∇s は表面勾配を意味し，∇s = P · ∇ である．

ρ = ρ1 φ + ρ2 (1 − φ),と評価できるので，ここでは，式 (6) より，(14)3.2 圧力駆動される赤血球の挙動
円柱内の流体を圧力駆動し，その中の赤血球形状をした弾性膜の挙動を解析する．ここで，弾性膜には Skalakモデルを用い，円柱はバルク弾性体 (Neo-Hookean) で与える．バルクの流体・構造体連成解析手法には，杉山ら(12,13,5,7)の方法を用い，本提案手法とカップリングする．
レイノルズ数はおよそ 1，赤血球内部には周囲流体の 5倍の粘性を与える．と与えられる．
オイラー型定式化では，界面は陰的に表現されている
ため，初期配置からの変形量を記述する量を予報変数として求める必要がある．Bs は，界面の法線ベクトルの情報が含まれているが，法線ベクトルは VOF 関数を用いて，∇φ,(15)n=−|∇φ|2Copyright c 2010 by JSFM

第 24 回数値流体力学シンポジウム
B5-5弾性膜の時間挙動を Fig.3 に示す．圧力勾配により駆動された流れにのり，弾性膜は変形しつつ移動しており，時間が進むにつれいわゆるパラシュート形状の様相をとってきている．本解析では，曲げ応力も考慮しているため，膜の尻尾部分の曲率変化が大きい部分が滑らかになっている．計算結果の最終ステップの主流方向に平均化した速度分布とポアズイユ分布を Fig.4 に示す．弾性膜が入っていることにより，流れ方向の抵抗が増大し速度分布が中央付近で一様になっている様子が確認できる．

(a) ca = 0.0125(b) ca = 0.05(a) t = 0(b) t = 0.05(c) t = 0.15(d) t = 0.25(e) t = 0.375(f) t = 0.5

Fig. 3: Time-evolutional behaviors of the biconcave
membrane capsules.

(c) ca = 0.2Fig. 1: Steady-state solution in shear ﬂow, with diﬀerent capillary numbers.

Fig. 4: The streamwisely-averaged velocity proﬁle forthe numerical ﬂow and Poiseuille ﬂow.
まとめと今後
本研究では，完全オイラー型定式化による流体・弾性膜連成解析手法を提案した．閉じた膜の内外を表現するために VOF 関数を導入し，更に，膜表面の応力を表現するために変形テンソルの時間発展を求めた．これにより，界面を陽に記述することなく膜の表面張力を計算することが可能となった．せん断場中の弾性膜の変形問題において本提案手法の妥当性を示し，また，圧力駆動される流れでの赤血球形状をした弾性膜の解析より，計算手法の強靭性を示した．今後は，より詳細な検証を行う

4.

Fig. 2: Time evolution of the deformation paramter
D(14) .

3

Copyright c 2010 by JSFM

第 24 回数値流体力学シンポジウム
B5-5とともに，より精密な赤血球および血流現象を再現するために，膜の透過現象，拡散現象を記述するモデルをオイラー型定式化で取り入れていく．
参考文献(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(1980) 251–267.

(2) C.D. Eggleton and A.S. Popel, Large deformationof red blood cell ghosts in a simple shear ﬂow, Phys.
Fluids, 10 (1998) 2182–2189.
(3) F.H. Harlow and J.E. Welch, Numerical calculation of time-dependent viscous incompressible ﬂowof ﬂuid with free surface, Phys. Fluids, 8 (1965)
1834–1845.
(4) C.W. Hirt, B.D. Nichols, Volume of ﬂuid (VOF)
methods for the dynamics of free boundaries, J.
Comput. Phys., 39 (1981) 201–225.
(5) S. Ii, K. Sugiyama, S. Takeuchi, S. Takagi andY. Matsumoto, An implicit full Eulerian method
for the ﬂuid-structure interaction problem, Int.
J. Numer. Meth. Fluids, in press (2010), DOI:
10.1002/ﬂd.2460.
(6) S. Ii and F. Xiao, A simpliﬁed interface capturing scheme using a coutinuous function, Proc. ofSNA+MC 2010, (2010).
(7) 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.
(8) C.S. Peskin, Flow patterns around heart valves: anumerical method, J. Comput. Phys., 10 (1972)
252–271.
(9) C. Pozrikidis, Finite deformation of liquid capsules
enclosed by elastic membranes in simple shear ﬂow,J. Fluid. Mech., 297 (1995) 123–152.
(10) 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.
(11) 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.
(12) 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.
(13) K. Sugiyama, S. Ii, S. Takeuchi, S. Takagi and Y.
Matsumoto, A full Eulerian ﬁnite diﬀerence approach for solving ﬂuid-structure coupling problem, J. Comput. Phys., in press (2010), DOI:
10.1016/j.jcp.2010.09.032.
(14) G.I. Taylor, The deformation of emulsions in deﬁnable ﬁelds of ﬂows, Proc. R. Soc. Lond. A, 146(1934) 501–523.

4Copyright c 2010 by JSFM